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Preface 



All practical systems contain nonlinear dynamics. Control system development for 
these systems has traditionally been based on linearized system dynamics in 
conjunction with linear control techniques. Sometimes it is possible to describe the 
operation of systems by a linear model around its operating points. Linearized system 
can provide approximate behavior of the system. But in analyzing the overall system 
behavior, the resulting system model is inadequate or inaccurate. Moreover, the 
stability of the system cannot be guaranteed. However, nonlinear control techniques 
take advantage of the given nonlinear dynamics to produce high-performance designs. 

Nonlinear Control Systems represent a new trend of investigation during the last few 
decades. There has been great excitement over the development of new mathematical 
techniques for the control of nonlinear systems. Methods for the analysis and design of 
nonlinear control systems have improved rapidly. A number of new approaches, ideas 
and results have emerged during this time. These developments have been motivated 
by comprehensive applications such as mechatronic, robotics, automotive and air-craft 
control systems. 

The book is organized into eleven chapters that include nonlinear design topics such 
as Feedback Linearization, Lyapunov Based Control, Adaptive Control, Optimal 
Control and Robust Control. All chapters discuss different applications that are 
basically independent of each other. The book will provide the reader with 
information on modern control techniques and results which cover a very wide 
application area. Each chapter attempts to demonstrate how one would apply these 
techniques to real-world systems through both simulations and experimental settings. 

Lastly, I would like to thank all the authors for their excellent contributions in different 
applications of Nonlinear Control Techniques. Despite the rapid advances in the field, 
I believe that the examples provided here allow us to look through some main 
research tendencies in the upcoming years. I hope the book will be a worthy 
contribution to the field of Nonlinear Control, and hopefully it will provide the 
readers with different points of view on this interesting branch of Control Engineering. 



Dr. Meral Alhnay 

Kocaeli University, 
Turkey 
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1. Introduction 

In nature, most of the systems are nonlinear. But, most of them are thought as linear and the 
control structures are realized with linear approach. Because, linear control methods are so 
strong to define the stability of the systems. However, linear control gives poor results in 
large operation range and the effects of hard nonlinearities cannot be derived from linear 
methods. Furthermore, designing linear controller, there must not be uncertainties on the 
parameters of system model because this causes performance degradation or instability. For 
that reasons, the nonlinear control are chosen. Nonlinear control methods also provide 
simplicity of the controller (Slotine & Li, 1991). 

There are lots of machine in industry. One of the basic one is dc machine. There are two 
kinds of dc machines which are brushless and brushed. Brushed type of dc machine needs 
more maintenance than the other type due to its brush and commutator. However, the 
control of brushless dc motor is more complicated. Whereas, the control of brushed dc 
machine is easier than all the other kind of machines. Furthermore, dc machines need to dc 
current. This dc current can be supplied by dc source or rectified ac source. Three - phase ac 
source can provide higher voltage than one phase ac source. When the rectified dc current is 
used, the dc machine can generate harmonic distortion and reactive power on grid side. 
Also for the speed control, the dc source must be variable. In this paper, dc machine is fed 
by three - phase voltage source pulse width modulation (PWM) rectifier. This kind of 
rectifiers compared to phase controlled rectifiers have lots of advantages such as lower line 
currents harmonics, sinusoidal line currents, controllable power factor and dc - link voltage. 
To make use of these advantages, the filters that are used for grid connection and the control 
algorithm must be chosen carefully. 

In the literature there are lots of control methods for both voltage source rectifier and dc 
machine. References (Ooi et al, 1987; Dixon&Ooi, 1988; Dixon, 1990; Wu et al, 1988, 1991) 
realize current control of L filtered PWM rectifier at three - phase system. Reference (Blasko 
& Kaura, 1997) derives mathematical model of Voltage Source Converter (VSC) in d-q and 
a-p frames and also controlled it in d-q frames, as in (Bose, 2002; Kazmierkowski et al, 
2002). Reference (Dai et al., 2001) realizes control of L filtered VSC with different decoupling 
structures. The design and control of LCL filtered VSC are carried out in d-q frames, as in 
(Lindgren, 1998; Liserre et al., 2005; Dannehl et al., 2007). Reference (Lee et al., 2000; Lee, 
2003) realize input-output nonlinear control of L filtered VSC, and also in reference 
(Komurciigil & Kiikrer, 1998) Lyapunov based controller is designed for VSC. The feedback 
linearization technique for LCL filtered VSC is also presented, as in (Kim & Lee, 2007; Sehirli 
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& Altinay, 2010). Reference (Holtz, 1994) compares the performance of pulse width 
modulation (PWM) techniques which are used for VSC. In (Krishnan, 2001) the control 
algorithms, theories and the structure of machines are described. The fuzzy and neural 
network controls are applied to dc machine, as in (Bates et al., 1993; Sousa & Bose, 1994). 

In this chapter, simulation of dc machine speed control which is fed by three - phase voltage 
source rectifier under input - output linearization nonlinear control, is realized. The speed 
control loop is combined with input-output linearization nonlinear control. By means of the 
simulation, power factor, line currents harmonic distortions and dc machine speed are 
presented. 

2. Main configuration of VSC 

In many industrial applications, it is desired that the rectifiers have the following features; 
high-unity power factor, low input current harmonic distortion, variable dc output voltage 
and occasionally, reversibility. Rectifiers with diodes and thyristors cannot meet most of 
these requirements. However, PWM rectifiers can provide these specifications in 
comparison with phase-controlled rectifiers that include diodes and thyristors. 

The power circuit of VSC topology shown in Fig.l is composed of six controlled switches 
and input L filters. Ac-side inputs are ideal three-phase symmetrical voltage source, which 
are filtered by inductor L and parasitic resistance R, then connected to three-phase rectifier 
consist of six insulated gate bipolar transistors (IGBTs) and diodes in reversed parallel. The 
output is composed of capacitance and resistance. 
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Fig. l.L filtered VSC 



3. Mathematical model of the VSC 

3.1 Model of the VSC in the three-phase reference frame 

Considering state variables on the circuit of Fig.l and applying Kirchhoff laws, model of 
VSC in the three-phase reference frame can be obtained, as in (Wu et al., 1988, 1991; Blasko 
& Kaura, 1997). 

The model of VSC is carried out under the following assumptions. 

• The power switches are ideal devices. 

• All circuit elements are LTI (Linear Time Invariant) 

• The input AC voltage is a balanced three-phase supply. 



Application of Input-Output Linearization 



For the three-phase voltage source rectifier, the phase duty cycles are defined as the duty 
cycle of the top switch in that phase, i.e., d a = d(Si), db= d(Ss), d c = d(Ss) with d representing 
duty cycle. 
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dt 
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dt 

di c 
dt 
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This model in equations (1) - (4) is nonlinear and time variant. Using Park Transformation, 
the ac-side quantities can be transformed into rotating d-q frame. Therefore, it is possible to 
obtain a time-invariant system model with a lower order. 

3.2 Coordinate transformation 

On the control of VSC, to make a transformation, there are three coordinates whose relations 
are shown by Fig 2, that are a-b-c, a-p and d-q. a-b-c is three phase coordinate, a-p is 
stationary coordinate and d-q is rotating coordinate which rotates a speed. 
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Fig. 2. Coordinates diagram of a-b-c, a-p and d-q 
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The d-q transformation is a transformation of coordinates from the three-phase stationary 
coordinate system to the d-q rotating coordinate system. A representation of a vector in any 
n-dimensional space is accomplished through the product of a transpose n-dimensional 
vector (base) of coordinate units and a vector representation of the vector, whose elements 
are corresponding projections on each coordinate axis, normalized by their unit values. In 
three phase (three dimensional) space, it looks like as in (5). 



X n hr — a n 



c u ] 



(5) 



Assuming a balanced three-phase system, a three-phase vector representation transforms to 
d-q vector representation (zero-axis component is 0) through the transformation matrix T, 
defined as in (6). 



T = 



2 2 

cos(cot) cos(cot — -it) cos(cot + -n) 

2 2 

— sin((x)t) — sin((x)t — n) — sm(cot + -n) 



(6) 



In (6), cd is the fundamental frequency of three-phase variables. The transformation from X abc 
(three-phase coordinates) to X d q(d-q rotating coordinates), called Park Transformation, is 
obtained through the multiplication of the vector X abc by the matrix T, as in (7). 



X dq — T.X a 



(7) 



The inverse transformation matrix (from d-q to a-b-c) is defined in (8). 



— s'm(a>t + -it) 



cos(cot) 
r _ t cos(<ut - \ii) -sm(cot-\n) 
cos(&)t + 2 -n) _ sin(&)t + l n -) 

The inverse transformation is calculated as in (9). 

%abc ~ T -Xdq 



(8) 



(9) 



3.3 Model of the VSC in the rotating frame 

Let x and u be the phase state variable vector and phase input vector in one phase of a 
balanced three-phase system with the state equation in one phase as in (10). 



X = Ax + Bu (10) 

Where A and B are identical for the three phases. Applying d-q transformation to the three- 
phase system, d-q subsystem with d and q variables is obtained (xd-x q and u<i-u q ). The 
system equation in (10) becomes as in (11) (Mao et al., 1998; Mihailovic, 1998). 



r A coh [ x di rB 0] [""di 
i-col A 1 UJ Lo B\ KJ 



(11) 



Where I is the identity matrix and is a zero matrix, both having the same dimension as x. 
(11) can transform any three-phase system into the d-q model directly. 
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When equations (1) - (4) are transformed into d-q coordinates, (12) - (14) are obtained, as in 
(Blasko & Kaura, 1997; Ye, 2000; Kazmierkowski et al, 2002). 
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dt 
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~~ yVdcdd ' 
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(12) 
di q R 1 U q 

■dr = -L t <- at *-L v « !d <-T (13) 

dV dc 11 MA . 

-j£- = -^{ldd d + Iqdq) - fjldc ( 14 ) 

Where i d and i q are the d-q transformation of i a , i b and i c . v d and v q are the d-q 
transformation of v a , v b and v c . d d and d q are the d-q transformation of d a , d b and d c . 

4. Input-output feedback linearization technique 

Feedback linearization can be used as a nonlinear design methodology. The basic idea is 
first to transform a nonlinear system into a (fully or partially) linear system, and then to use 
the well-known and powerful linear design techniques to complete the control design. It is 
completely different from conventional linearization. In feedback linearization, instead of 
linear approximations of the dynamics, the process is carried out by exact state 
transformation and feedback. Besides, it is thought that the original system is transformed 
into an equivalent simpler form. Furthermore, there are two feedback linearization methods 
that are input-state and input-output feedback linearization (Slotine & Li, 1991; Isidori, 1995; 
Khalil, 2000; Lee et al., 2000; Lee, 2003). 

The input-output feedback linearization technique is summarized by three rules; 

• Deriving output until input appears 

• Choosing a new control variable which provides to reduce the tracking error and to 
eliminate the nonlinearity 

• Studying stability of the internal dynamics which are the part of system dynamics 
cannot be observed in input-output linearization (Slotine & Li, 1991) 

If it is considered an input-output system, as in (15)-(16); 

X - /O) + g(x)u (15) 

V = Hx) (16) 

To obtain input-output linearization of this system, the outputs y must be differentiated 
until inputs u appears. By differentiating (16), equation (17) is obtained. 

dh 

y = fa 1/00 + 500"] = L f h(x) + L g h(x)u (17) 

In (17), Lfh and L g h are the Lie derivatives of f(x) and h(x), respectively and identified in (18). 



dh dh 

-f{x), L g h(x)=- i 



L / ft O0 = — fix) , L g h(x) = —g(x) (18) 
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If the k is taken as a constant value; k. order derivatives of h(x) and 0. order derivative of 
h(x) are shown in (19) - (20), respectively. 



L k f h(x) = Lf]}f x h{x) 



d{fif x K) 



u f 



■> dx 

V}h(x) - h(x) 



fix) 



(19) 



(20) 



After first derivation, If L g h is equal to "0", the output equation becomes y — Lfh(x). 
However, it is independent from u input. Therefore, it is required to take a derivative of 
output again. Second derivation of output can be written in (23), with the help of (21)-(22). 



dL f h 



d(L f h) 
LglfKx) = gx g(x) 

L 2 f h(x) = L f L f h{x) = \[_ f(x) 



dx 



'f ? 

[fix) + gix)u\ — L,h(x) + L g Lfh(x)u 



(21) 
(22) 
(23) 



dx 

If L g Lfh(x) is again equal to "0", y is equal to L?/i(x)and it is also independent from u input 
and it is continued to take the derivation of output. After r times derivation, if the condition 
of (24) is provided, input appears in output and (25) is obtained. 



L g .L r f ' %(x) * 



yp^L^+^i^Ly- 1 ^)^ 



(24) 
(25) 



Applying (25) for all n outputs, (26) is derived. 
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(26) 



E(x) in (27) is a decoupling matrix, if it is invertible and new control variable is chosen, 
feedback transformation is obtained, as in (28). 
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(28) 



Equation (29) shows the relation between the new inputs v and the outputs y. The input- 
output relation is decoupled and linear (Lee et al., 2000). 
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If the closed loop error dynamics is considered, as in (30) - (31), (32) defines new inputs for 
tracking control. 



eX + k 



l(r-2) e l 



■ + ■■• + tiie, 1 + /c ln ei 



p r A- h p r ~^ 

c n "■" ls -n(r-l) c n 



+ ■■■ + k 21 el + k 20 e n 







(30) 



vy - y 

r-l k . .,,1 . 



-ki(r-i)y r kn(r-i)y - fcio(yi - yl) 

~k n (r-i)y r — k 2 i(r-i)y — k 20 (y n —yn) 



(31) 



(32) 



k values in equations show the constant values for stability of systems and tracking of y 
references, as in (Lee, 2003). 

5. The application of an input-output feedback linearization to the VSC 

The state feedback transformation allows the linear and independent control of the d and q 
components of the line currents in VSC by means of the new inputs u d and u q . 

For unity power factor, in (12 - 14) u d — V m and u q — are taken, so mathematical model of 
this system is derived with (33-35), as in (Komiirciigil & Kiikrer, 1998; Lee, 2003). 



diii 
dt 
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— i d + a>i q --V dc d d — — 
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— -- I l q -coi d - I V dc d q 
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If (33-35) are written with the form of (15), (36) is derived. 
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(36) 



The main purpose of this control method is to regulate V dc voltage by setting i d current and 
to provide unity power factor by controlling i q current. Therefore, variables of y outputs 
and reference values are chosen as in (37). 
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Differentiating outputs of (37), (38) is obtained. The order of derivation process, finding a 
relation between y outputs and u inputs, is called as relative degree. It is also seen that the 
relative degree of the system is '1'. 
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Fig. 3. Input-output feedback linearization control algorithm of VSC 
When (38) is ordered like (28), (39) is obtained. 
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After taking inverse of matrix (39) and adding new control inputs from (40), (41) is obtained. 
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Control algorithm is seen in Fig.3. For both L and LCL filtered VSC, the same control 
algorithm can be used. Providing the unity power factor, angle values are obtained from 
line voltages. This angle values are used in transformation of a-b-c to d-q frames. Line 
currents which are transformed into d-q frame, are compared with d-q reference currents, d 
axis reference current i^ref i s obtained by comparison of dc reference voltage V re f and actual 
dc voltage V dc . On the other hand, q axis reference current i qre f is set to '0' to provide unity 
power factor. And by using (41), switching functions of d-q components are found. Then, 
this d-q switching functions are transformed into a-b-c frames and they are sent to PWM 
block to produce pulses for power switches. 

Providing the control of internal dynamics, in dc controller square of V re f and V dc are used 
(Lee, 2003). 



6. DC machine and armature circuit model 

Electrical machines are used for the conversion of electric power to mechanical power or 
vice versa. In industry, there are wide range of electrical machines that are dc machines, 
induction machines, synchronous machines, variable reluctance machines and stepping 
motors. The Dc machines can be classified as a brushless and brushed dc machines. 
Furthermore, the advantage of brushed dc machines is the simplicity with regard to speed 
control in the whole machines. However, the main disadvantage of this kind of machines is 
the need of maintenance because of its brushes and commutators. 

Fig. 4 shows the basic structure of brushed dc machines. Basic components of dc machines 
are field poles, armature, commutator and brushes (Fitzgerald et al., 2003). 




Fig. 4. Dc machine 

Field poles produce the main magnetic flux inside of the machines with the help of the field 
coils which are wound around the poles and carry the field current. Some of the dc 
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machines, the magnetic flux is provided by the permanent magnet instead of the field coils. 
In Fig. 5, the field coils and field poles of dc machines are shown (Bal, 2008). 




Fig. 5. Field coils and field poles of a dc machine 

The rotating part of the dc machine is called as an armature. The armature consists of iron 
core, conductors and commutator. Besides, there is a shaft inside of armature that rotates 
between the field poles. The other part of the machine is commutator which is made up of 
copper segments and it is mounted on the shaft. Furthermore, the armature conductors are 
connected on the commutator. Another component of dc machine is brushes. The brushes 
provide the electric current required by armature conductors. In dc machine to ensure the 
rotation of the shaft, the armature conductors must be energized. This task is achieved by 
brushes that contact copper segments of commutator. Also, the brushes generally consist of 
carbon due to its good characteristic of electrical permeability. Fig. 6 shows the armature, 
commutator and the brushes (Fitzgerald et al., 2003; Bal, 2008). 




Fig. 6. a) armature, b) commutator and c) brushes of a dc machine 
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To produce the main flux, the field must be excited. For this task, there are four methods 
which are separately, shunt, series and compound to excitation of dc machines and are 
shown in Fig. 7. However, separately excited dc machine is the most useful method because 
it provides independent control of field and armature currents. Therefore, this structure is 
used in this chapter (Krishnan, 2001; Fitzgerald et al., 2003). 



la La Ra 
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b) 



Series^ 
field i 1 ^ 

Ra,L a 



V 



IfV 



Shunt 
field f 

IfRf 



c) I 

Fig. 7. Excitation methods of dc machine a) separately, b) shunt, c) series, d) compound 
excitation 

There are two basic speed control structure of dc machine which are armature and field, as 
in (Krishnan, 2001). The armature circuit model of dc machine is shown in Fig. 8. 





Fig. 8. Armature circuit model of dc machine 
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The mathematical model of armature circuit can be written by (42). 

v = e + R a I a +L a ^- (42) 

In steady state, — - part is zero because of the armature current is constant. The armature 
model is then obtained by (43), (Krishnan, 2001). 

v = e + R a I a (43) 

e = K<P f a> m (44) 

(43) is written in (44), (45) is derived. 

a> m ~ (45) 

l f 

The speed of dc machine depends on armature voltage and field current, as shown in (45). In 
field control, the armature voltage is kept constant and the field current is set. The relation 
between speed and field current is indirect proportion. However, in armature control, the 
relation between armature voltage and speed is directly proportional. Furthermore, in 
armature control, the field current is kept constant and the armature voltage is set. 

In this chapter, the armature control of dc machine is realized. 

The speed control loop is added to nonlinear control loop. Firstly, the actual speed is 
compared with reference speed then the speed error is regulated by PI controller and after 
that its subtraction from armature current, the reference current is obtained. The reference 
current obtained by speed loop, is added to nonlinear control loop instead of reference i d 
current, which is obtained by the comparison of the square of reference voltage and actual 
voltage. 



In Fig. 9, speed control loop is shown. 



W ref+ ~A e 



i+rV 



dref 



4 i 

Fig. 9. Armature speed control loop of dc machine 

7. Simulations 

Simulations are realized with Matlab/Simulink. Line voltage is taken 220 V, 60 Hz. The 
switching frequency is also chosen 9 kHz. L filter and controllers parameters are shown in 
Table.l and Table.2. 

Simulation diagram is shown in Fig.10. By simulation, steady-state error and settling time of 
dc motor speed, harmonic distortions and shapes of line currents and unity power factor are 
examined. 
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Passive Components 


L Filter 


De-Link 


L(H) 


R(Q) 


Cdc (P F ) 


0.0045 


5.5 


2200 



Table 1. Values of L filter components 



Controllers 


Speed Controller 


Input-output current controller 


*v 


Ki 


k (10 3 ) 


K (10 s ) 


10 


0.01 


30 


50 



Table 2. Values of controllers 
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Fig. 10. Simulation diagram of dc machine controller in Simulink 
Fig. 11 shows the structure of input-output controller diagram. 
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Fig. 11. Input-output controller diagram 
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Equation (41) is written in the block of input-output controller which is shown in Fig. 12. 
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Fig. 12. Input - output controller 

Fig. 13 shows the speed controller of dc machine. 
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Fig. 13. Speed controller of dc machine 

Fig.14 shows the dc machine speed. Reference speed value is changed from 150 rad/s to 
200rad/s at 0.5 s. Settling time to the first reference is shorter than 0.15 s, but settling time of 
second reference is 0.1 s. 

Fig. 15 shows the steady-state error of dc machine speed. It is seen that the steady - state 
error changes between ±2 rad/s. 

The one phase voltage and current is shown in Fig. 16. It is also seen that unity power factor 
is obtained but not as desired. 



Fig. 17 shows the line currents. The shapes of line currents are sinusoidal. 
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Fig. 16. One phase voltage and current 
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Fig. 18 shows the harmonic distortions of line currents. Line currents include high order 
harmonic contents. However, total harmonic distortion value (THD) is under the value that 
is defined by standards. THD of line currents are %1.34, %1.71 and %2.84. 

Fundamental (60Hz) = 35.63 , THD= 1.34% 
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Fig. 18. Harmonic distortions of line currents 
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8. Conclusion 

In this chapter, simulation of dc machine armature speed control is realized. Dc machine is 
fed by voltage source rectifier which is controlled input - output linearization nonlinear 
control method. Furthermore, for the speed control, dc link voltage is regulated by the dc 
machine speed control loop. The control algorithm of voltage source rectifier and dc motor 
speed are combined. The required reference I d current for voltage source rectifier is obtained 
by speed control loop. Simulations are carried through Matlab/Simulink. By means of the 
simulation results, the speed of dc machine, line currents harmonic distortions and power 
factor of grid are shown. It is shown that the voltage source rectifier with dc machine as a 
load provides lower harmonic distortion and higher power factor. Furthermore, dc machine 
speed can be regulated. 
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1. Introduction 

Many power electronic system designs focus on energy conversion circuit parameters rather 
than controller parameters which drive the circuits. Controllers must be designed on the 
basis of circuit models, which are generally nonlinear systems (Brockett & Wood (1974)), in 
order to improve performance of controlled systems. The performance of controlled systems 
depend on nominal models to design the controllers. More broad class of models controllers 
are designed for, better control performance may be given. Many works (e.g. Kassakian 
et al. (1991)) design controllers for linearized models because it is not easy to concretely 
design controllers for the nonlinear models. Controller design in consideration of nonlinear 
models has been discussed since a work by Banerjee & Verghese (2001) because a research on 
nonlinear controller design has grown in those times. 

This chapter systematically designs a robust controlled power converter system on the 
basis of its nonlinear model. Concretely, a two-stage power factor correction converter, 
that is a forward converter (FC) with power factor corrector (PFC), is designed. The 
systematic controller design clearly analyzes the behavior of nonlinear system to improve 
the performance. 

A work by Orabi & Ninomiya (2003) analyzes a stability of single-stage PFC for variations 
of controller gain on the basis of its nonlinear model. On the basis of the work by Orabi & 
Ninomiya (2003) that regards a load of PFC as constant, a work by Dranga et al. (2005) for a 
two-stage PFC focuses on a point that a load of PFC part is not pure resistive and analyzes a 
stability of the converter. 

A work by Sasaki (2002) discusses an influence between a FC part and a PFC part in 
a two-stage PFC. A work by Sasaki (2009) clearly shows that a source current reference 
generator plays an important role in a synthesis for a single-stage PFC. For the works by 
Sasaki (2002; 2009), this chapter shows how to decide synthesis parameters of robust linear 
and nonlinear controllers for a two-stage PFC in more detail. 

The controller synthesis step, that is this chapter, is organized as follows. First, the converter is 
divided into two parts which consists of a FC part and a PFC part by considering an equivalent 
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circuit of transformer. The two parts depend on each other and are nonlinear systems. The 
FC part has an apparent input voltage which depends on an output voltage in the PFC part. 
On the other hand, the PFC part has an apparent load resistance which depends on an input 
current in the FC part. Second, the two parts of converter are treated as two independent 
converters by analyzing steady state in the converter and deciding a set point. Then, the 
above input voltage and load resistance are fixed on the set point. Controllers are designed 
for the two independent converters respectively. For the FC part which is a linear system at 
the second stage of converter, a robust linear controller is designed against variations of the 
apparent input voltage and a load resistance. For the PFC part which is a bilinear system, 
that is a class of nonlinear system, at the first stage, a robust nonlinear controller is designed 
against variations of the apparent load resistance that mean dynamic variations of the FC part 
at the second stage. Finally, computer simulations demonstrate efficiencies of the approach. 
It is also clarified that consideration of nominal load resistance for each part characterizes a 
performance of the designed robust controlled system. 



2. Two-stage power factor correction converter 

A controlled system for a two-stage power factor correction converter, that is a forward 
converter (FC) with power factor corrector (PFC), is systematically constructed as shown in 
Fig.l. The systematic controller design clearly analyzes the behavior of controlled converter 
system. In this section, first, an equivalent circuit of transformer divides the converter into two 
parts with FC and PFC. Averaged models for the two parts are derived respectively, which 
depend on each other. Second, steady state in the two parts is analyzed, which is used as a set 
point. Finally, each of parts is treated as an independent converter respectively. 



''2 h^ L 2 
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Fig. 1. Two-stage converter ; forward converter (FC) with power factor corrector (PFC) 



2.1 Nonlinear averaged models 

An equivalent circuit of transformer derives a circuit as shown in Fig. 2 from Fig. 1. A FC part 

has a variable dc source voltage which depends on a switch S2. A PFC part has a variable load 

fc 
resistance which depends on a switch Sj. Fig. 2 gives two nonlinear averaged models (E-ca) 

and (£c,i ) for the FC and the PFC parts respectively, which depend on each other. 
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The FC model v^c/J is given as 



d 


'v{ 




1 

RCi 


1 


dt 


M_ 




1 

IT 


'1 
IT 





*>1 




f 









+ 


[07 


N 




_'!_ 




I 


LnJ 



n 



(i) 



Lyapunov-Based Robust and Nonlinear Control 
for Two-Stage Power Factor Correction Converter 



23 




Fig. 2. Two converters ; FC and PFC 
and the PFC model (E^ c ) is 
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where V\ and v 2 are averaged capacitor voltages, i\,i 2 averaged inductor currents, C\,C 2 
capacitances, L\,L 2 inductances, r\,r 2 internal resistances, R a load resistance, N a turns ratio, 
v s a source voltage, fi\,}i 2 averaged switching functions given by < yl\ < 1 and < fl 2 < 1. 

2.2 DC components of steady state 

DC components of steady state in two parts of converter are derived, which are used as a 
set point for controller design. Given an average full-wave rectified voltage V s := (co/ln) 

— — i— 

Jr." \v s (t)\dt = 2\/2V e /n of a source voltage v s = \/2V e smcvt, and specify averaged 

switching functions y.\ = fi\ s , fl 2 = fl 2s which are called a set point, then dc components 

of voltages and currents of steady state are given by 
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The above steady states in two parts depend on each other through (7). It is easily clarified by 
assuming r\ = r 2 = that the equations (3), (4) give a steady state in buck converter and the 
equations (5), (6) give a steady state in boost converter. 

2.3 Two independent models 

For two parts of converter with FC and PFC, the steady state analysis derives two independent 
converters. It means that the FC has an apparent source voltage £(f'2s) an d the PFC has an 
apparent load resistance -R2(Fi s ) which are fixed on a set point respectively. 

Then, two independent averaged models of the converters are given as 

fc 
a linear FC model C^aq) of the form 
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and a nonlinear PFC model (£//, ) of the form 
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The FC is modeled on a linear system and the PFC is modeled on a bilinear system that is a 
class of nonlinear system (Brockett & Wood (1974); Mohler (1991)). 

Here derived is an amplitude of source current of steady state in the PFC, which is used for 
source current reference generator in Section 3.3. Given a source voltage v s = \/2V e sin cvt [V], 
assume an output voltage Vi = Vi r [V], then dc components of steady state in the PFC gives a 
source current ii = \fll e sin wf [A], which is sinusoidal and in phase with the source voltage, 
given by 

(10) 



le 



In 



-2r 



4, i 



riRiiVis) 



This decision process is similar to that is discussed in works by Escobar et al. (1999); Escobar 
et al. (2001). 



3. Robust controller design 

Controllers for two parts with FC and PFC are designed on the basis of two independent 
converter models (£^ ) and (£^g ) respectively. First, the models (S^m) and (£^g ) are moved 
to set points given in Section 2.2. Next, controller design models, that are called as generalized 
plants, are derived which incorporate weighting functions in consideration of controller 
design specifications. Finally, feedback gains are given to guarantee closed-loop stability and 
reference tracking performance. 
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3.1 Models around set points 

Converter models around set points are derived from the averaged models (£%) an d (^ao )• 
Moving the state to a specified set point given by V\ := V\ — V\ s , i\ := l\ — i\ s , fa := }i\ — }i\ s , 
V2 '■= &2 — Vjs, 12 '■= h ~ hs, fa '■= fa ~ fas derives averaged models around the set point 
respectively for the two parts, which are given by 

a FC model (E\, ) of the form 
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and a PFC model (L. ) of the form 
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fi2s- The following discussion constructs 



controllers on the basis of the averaged models (E-\ ) and (Ej7 ) around the set point. 



3.2 Linear controller design for FC 

A robust controller for the FC is given by linear H 00 control technique because the averaged 

fc 
model (L' A ) is a linear system. The control technique incorporates weighting functions 

W B i(s) := k v i/(s + £i) to keep an output voltage constant against variations of load resistance 

and apparent input voltage as shown in Fig. 3, where k v \, £\ are synthesis parameters. Then, 

the linear H°° control technique gives a linear gain Ki of the form 
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and xi, denotes state 



tp * w j , ->2 ■ [Op U J ' '-'12 
of weighting function W v i(s) as shown in Fig.3. In the figure, matrices W e j and W u i are 
weighting coefficients of a performance index in the linear H°° control technique and used 
only in the controller synthesis. 

The matrix Yf c that constructs the gain K-y is given as a positive-definite solution satisfying a 
Lyapunov-based inequality condition 
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Fig. 3. Block diagram of linear controlled system for FC 
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c{ c yf c -i 

where the matrix AJ C ,B^ , C[ are coefficients of a generalized plant given by Fig.3 and ■yf is 
a synthesis parameter in the linear H°° control technique. 



3.3 Nonlinear controller design for PFC 

A robust nonlinear controller design for the PFC, that is a nonlinear system, consists of the 
following two steps. First, a nonlinear gain is given to guarantee closed loop system stability 
and reference tracking performance for source current and output voltage against variations of 
the apparent load resistance. Second, a source current reference generator is derived to adjust 
an amplitude of current reference to variations of source voltage and load resistance. At the 
first step, as shown in Fig. 4 incorporated are a weighting function W{2.(s) := kQ(v%/(s + 
2£,Wj2S + WJ2) for a source current to be sinusoidal and be in phase with a source voltage and 
a function W 2(s) := k V 2/(s + £2) for an output voltage to be kept constant against those 
variations, where kj2, Z,, U3ii, k V 2, £2 are synthesis parameters. Then, a nonlinear H 00 control 
technique in a work by Sasaki & Uchida (1998) gives a nonlinear gain of the form 
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as shown in Fig.4, where xPf c = \ x Pf cT x P/ cT ] , D[{ c = [0 T W,f 2 ] T , B p 2 fc (xPf c ) 
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' B p2i = [° rj] ' B p22 C = [ - tS °] and x w C denotes state of 

weighting functions. In the figure, matrices W e 2 and W„2 are weighting coefficients of a 
performance index in the nonlinear H°° control technique. A block named as Generator is 
not used for the synthesis and is discussed at the following second step. 
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Fig. 4. Block diagram of nonlinear controlled system for PFC 

The matrix yP/ c is given as a positive-definite solution satisfying a state-depended 
Lyapunov-based inequality condition 
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where x is a state of a generalized plant given by Fig.4, matrices AP' C ,B? ,5? ( X )/Cf /Dfo 
are coefficients of the plant and 7^ is a synthesis parameter in the nonlinear H°° control 
technique. 

For any state x, which is current and voltage in a specified domain, the matrix YV> C satisfying 
the inequality (18) is concretely given by solving linear matrix inequalities at vertices of a 
convex hull enclosing the domain as shown in a work by Sasaki & Uchida (1998). 

Next, given is a mechanism to generate a source current reference iir- As shown in Fig.4 
the reference generator consists of a feedforward loop given by steady state analysis in 
Section 2.3 and a feedback loop with a voltage error amplifier. The amplifier K(s) is given 
by -K(s) = kp + kj/s + kps where kp,kj and kp are constant parameters decided by system 
designers as discussed in a work by Sasaki (2009). The feedback loop is the same structure 
as a conventional loop used in many works (e.g., Redl (1994)). Note that the amplifier K(s) 
works only for variations from the nominal values in the circuit, because the feedforward loop 
gives the effective value of the source current as shown in Section 4.4. 



4. Computer simulations 

This section finally shows efficiencies of the approach through computer simulations. It 
is also clarified that consideration of nominal load resistance for each part characterizes a 
performance of the designed controlled system. A software package that consists of MATLAB, 
Simulink and LMI Control Toolbox is used for the simulations. 
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Parameters of the circuit shown in Fig.l are given by Table 1. 



r-i 1 [mfl] 



Lj 10 [mH] 



Ci 450 [f(F] 



N 1/36 



r 2 1 [mfl] 



L 2 1 [mH] 



C 2 1000 [jiF] 



Table 1. Circuit parameters of two-stage power factor correction converter as shown in Fig.l 

4.1 Design specification 

Control system design specification for the two-stage power factor correction converter is 
given by 

(1) An effective value of source voltage is 82 ~ 255 volts, and its frequency is 45 ~ 65 hertz ; 

(2) An output voltage is kept be 5 volts, and its error of steady state is within ±0.5%; 

(3) An output current is ~ 30 amperes ; 

(4) A source current is approximately sinusoidal and is phase with source voltage. 

The above specifications are treated for controller design as the following respective 
considerations ; 

(1) A nonlinear gain for PFC is designed for a source voltage whose nominal effective value 
is V e = (82 + 255) /2 = 168.5 volts (then, average full-wave rectified voltage is V s = 151.7 
volts). Moreover, the gain is designed to be robust against source voltage variations by 
treating that as a disturbance wJ as shown in Fig.4. 

(2) Incorporated is a weighting function W„i(s) whose magnitude is high at low frequencies. 

(3) A load resistances R varies between 1/6 ~ oo ohms at an output voltage of 5 volts. 

(4) Incorporated is a weighting function Wq(s) whose magnitude is high at frequencies in 
the range of 45 ~ 65 hertz for the source current to be approximately sinusoidal at the 
frequencies. 

4.2 Nominal load resistance 

Now, a nominal value of load resistance R need be decided to design controllers. The nominal 
value influences a performance of closed-loop system controlled by the designed controllers. 
Root loci of linearized converters for the two parts as shown in Figs. 5 and 6 , which are 
pointwise eigenvalues of the system matrices for variations of load resistance R, show that the 
larger the load resistance R is, the more oscillatory the behavior of FC is and the slower the 
transient response in PFC is. Controllers, here, are constructed in consideration of undesirable 
behavior of the converter system. Therefore, a FC controller is designed for a system with 
oscillatory behavior, that is for a large load resistance, so that an output voltage does not 
oscillate. A PFC controller is designed for a system with fast response, that is for a small load 
resistance, so that a source current is not distorted by the controller responding well to source 
voltage variations. 
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Fig. 5. Root locus of coefficient A„ of linear model for FC 
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Fig. 6. Root locus of coefficient A?J of linearized model for PFC 
4.3 Linear controller design parameters for FC 

A FC controller is designed for a large nominal load resistance R = 100 ohms. Consider a set 
point \i\ s = 0.5 for a capacitor voltage of 02s = 360 volts. DC components of steady state are 
given by V\ s = 4.99 volts and i\ s = 0.0499 amperes. For a model around the above set point, 
incorporated are weighting functions with parameters given by K v i = 20, £\ = 0.001 W e i = 5, 
W„i = 1.5 and yJ c = 0.95 where jf c denotes control performance level given in the linear H°° 
control technique. The weighting function W v \(s) with the above parameter gives bode plots 
whose characteristics is like an integrator as shown in Fig. 7. 

Then, the matrix Y> c that gives a linear gain (15) is obtained as 



yfc 



2.00 x 10 4 -1.55 x 10 3 5.05 x 10 2 ' 
-1.55 x 10 3 7.68 x 10 2 2.59 x 10 1 
5.05 x 10 2 2.59 x 10 1 2.51 x 10 1 



(19) 
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Fig. 7. Bode plots of weighting function W v i(s) for output voltage tracking in FC controller 
design 

4.4 Nonlinear controller design parameters for PFC 

A PFC controller is designed for a small nominal load resistance R = 10 ohms. Consider a 
set point p2s = { v 2r ~ V s )/V2 r = 0.579 to make a rectified voltage »2r = 360 volts. Then, an 
apparent load resistance in the PFC is given by 7?2(^ls) = 51845.2 ohms for a set point fii s = 
0.5 in the FC. Therefore, dc components of steady state are given by »2s = 359.9 volts and i^ = 




Fig. 8. Bode plots of weighting function Wq(s) for achieving unity power factor in PFC 
controller design 
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Fig. 9. Bode plots of weighting function W v2 (s) for capacitor voltage regulation in PFC 
controller design 

0.0165 amperes. For a model around the above set point, incorporated are weighting functions 
with parameters given by £ = 0.001, w i2 = 1307T, K i2 = 300 /cj a , K v2 = 0.002, e 2 = 0.001, 
W e2 = diag [10~ 7 , 1], W u2 = 20 and yPf c = 0.98 where yPJ c denotes control performance level 
given in the nonlinear H°° control technique. The weighting function Wj 2 (s) with the above 
parameter gives bode plots that is weighted at frequencies given by the design specification 
(4). The function W v2 (s) has low gains at all frequencies because a voltage tracking fed-back 
by a controller is not important for the PFC design. It also leads that the W e2 is set such that 
the tracking performance for source current is more weighted than that for output voltage. 

Then the matrix YPf c that gives a nonlinear gain (17) is obtained as 

2.77 x 10 4 -5.56 x 10 4 2.93 x 10° 4.70 x 10 3 -1.72 x 10 4 
-5.56 x 10 4 1.89 x 10 5 8.86 x 10" 1 -5.57 x 10 3 4.03 x 10 4 
2.93 x 10° 8.86 x 10" 1 



ypf c 



(20) 



9.71 x 10 8 -2.60 x 10~ 2 1.89 x 10" 1 
4.70 x 10 3 -5.57 x 10 3 -2.60 x 10~ 2 1.10 x 10 3 -2.98 x 10 3 
-1.72 x 10 4 4.03 x 10 4 1.89 x 10" 1 -2.98 x 10 3 1.24 x 10 4 

by considering a voltage variation of ±10 volts and a current variation of ±3 amperes around 
the set point (i.e., by considering 349.9 < v 2 < 369.9 and -2.9835 < i 2 < 3.0165) and solving 
a convex programming problem as shown in the work by Sasaki & Uchida (1998). 

Next, a source current reference generator is constructed with a voltage error amplifier given 
by kp = 8 x 10~ 2 , kj = 4, fcp = 1 x 10~ 5 as shown in Fig. 10. The gain is chosen in order to be 
low around twice the input line frequency. 
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Fig. 10. Bode plots of voltage error amplifier K(s) in PFC controller 



4.5 Simulation results 

Figs. 11-13 show behaviors of the nonlinear averaged models of FC (E M ) and PFC (£>c' A ) in 
the following cases ; 

(CIA load resistance R changes from 1000 to 0.25 ohms in steady state for a source voltage V s 
= \fl 100 sinl007rf volts ; 

(C2An efficient value V e of source voltage with 50 hertz changes from 100 to 85 volts for 2 
seconds in steady state with a load resistance R = 0.25 ohms. 

Figs. 11-13 show that the controllers works very well. 

Fig.12 shows behaviors by a different controller from that in Figs. 11 and 13, which is designed 
for desirable load resistances with a small nominal resistance R = 10 ohms for the FC and a 
large R = 100 ohms for the PFC. 

Therefore, in the FC the output voltage V\ in Fig. 12 oscillates more than in Fig. 11 at the large 
resistance R = 1000. Also, in the PFC the capacitor voltage Vi in Fig. 12 drops larger than 
in Fig. 11. Then the source current ij in Fig.12 raises faster than in Fig. 11 and then is more 
distorted. 

In the case as shown in Fig. 12, the above control techniques give the following matrices as 

2.92 x 10 4 -1.28 x 10 3 6.26 x 10 2 ~ 
Y? c = -1.28 x 10 3 8.07 x 10 2 3.57 x 10 1 (21) 

6.26 x 10 2 3.57 x 10 1 2.27 x 10 1 
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(a) output voltage V\ 



(b) output current i\ 



(c) FC controller output U\ 




(d) capacitor voltage v 2 



! '• 



time /[s] 

(e) source current 12 





(f) PFC controller output 112 



Fig. 11. (CI) Behaviors when a load resistance R changes from 1000 to 0.25 ohms in steady 
state 



and 



ypfc 



2.76 x 10 4 -5.56 x 10 4 2.93 x 10° 4.70 x 10 3 -1.72 x 10 4 

-5.56 x 10 4 1.89 x 10 5 8.85 x 10" 1 -5.56 x 10 3 4.03 x 10 4 

2.93 x 10° 8.85 x 10" 1 9.71 x 10 8 -2.60 x 10~ 2 1.89 x 10" 1 

4.70 x 10 3 -5.56 x 10 3 -2.60 x 10" 2 1.10 x 10 3 -2.98 x 10 3 

-1.72 x 10 4 4.03 x 10 4 1.89 x 10- 1 -2.98 x 10 3 1.24 x 10 4 



(22) 
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(a) output voltage V\ 



(b) output current i\ 





(c) FC controller output U\ 



(d) capacitor voltage v 2 





(e) source current 12 



(f) PFC controller output 112 



Fig. 12. (CI) Behaviors when a load resistance R changes from 1000 to 0.25 ohms in steady 
state where a controller is designed for a different nominal load resistance from that in Fig. 11 

Entries of the matrix Y> c in (21) are clearly different from those in (19). On the other hand, 
difference of the YPj c between the matrix (22) and (20) is much small. Those differences give 
the different behaviors as mentioned above. 

Fig. 12 demonstrates the discussion in Section 4.2 that controllers needs be designed in 
consideration of undesirable behavior of system. For a large load resistance the output voltage 
of the FC oscillates and for small load resistance the source current and the capacitor voltage 
of the PFC respond to variations larger than those in Fig. 11. 
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(a) output voltage V\ 



(b) output current i\ 




wwwwyviAfl/ 



(c) FC controller output U\ 



(d) capacitor voltage v 2 





(e) source current 12 



(f) PFC controller output 112 



Fig. 13. (C2) Behaviors when an efficient value of source voltage V e changes from 100 to 85 
volts for 2 seconds in steady state 

5. Conclusion 

For a two-stage power factor correction converter, a nonlinear controlled converter system 
was systematically designed on the basis of its nonlinear model. The systematic controller 
design clearly analyzed the behavior of nonlinear system to improve the performance. The 
synthesis step were clearly shown. It was clarified that a nominal load resistance characterizes 
the controller performance. Finally, computer simulations demonstrated efficiencies of 
the approach. The nonlinear controller synthesis was shown as a natural extension of a 
well-known Lyapunov-based linear controller synthesis. 
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1. Introduction 

Fluid systems where a liquid phase is dispersed in other liquid, as emulsions, are present in 
many industrial processes, technological applications and in natural systems. The flow of 
these substances shows a rheological behavior that depends on the viscosities ratio, the 
surface tension, surfactants, flow-type parameter and the coupled effects of the fluid 
structure and the kinematics properties of the flow, mainly in the non linear regimen, which 
is the reason because the dynamics of the fluid particles is an area of current research. As an 
approach to understand the physical properties of these fluids, studies on the deformation, 
break-up and coalescence of drops have been performed since the pioneering work of 
Taylor (1932). 

The deformation and breakup of liquid drops is dependent on the kinematics properties of 
the imposed flow, in particular of the second invariant of the rate of deformation tensor II2D 
and the flow-type parameter, a (see Astarita, 1979). In experimental studies of the dynamics, 
break-up and coalescence of drops, then it is necessary to be able to modify on demand the 
external flow field parameters causing the drop deformation. For example, Taylor (1934) use 
two flow devices, a Parallel Band Apparatus (PBA) and a Four-Roll Mill (FRM), which 
manually manipulated each of the rollers speeds to position a drop. Once the drop was in 
place, Taylor was able to track changes that occur on the drop shape as a function of the 
imposed flow, although for a short time, until the drop was ejected from its unstable 
position. 

Four-Roll Mills and Parallel Band cover a wide interval in the flow-type parameter; 
however, there is a gap between them that the flow fields generated by co-rotating Two-Roll 
Mill (TRM) geometries fill. PBA works for simple shear flow, corresponding to a flow-type 
parameter a = 0; FRM setup is more effective in the interval 0.4 < a < 1 [Yang et al. (2001)], 
whereas TRM is effective in the interval 0.03 < a < 0.3 [Reyes and Geffroy (2000b)]. 

In the case of the Four-Roll Mill, Bentley and Leal (1986a) have shown how to control -for 
long times- the position of drops within the flow field generated by a FRM by adjusting with 
a computer and in real-time the speed of rotation of each cylinder, and maintaining 
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simultaneously constant values for the elongational flow field parameters. This is achieved 
by modifying the angular velocity of the cylinders, which displaces the stagnation point to a 
new position, in order to drive the motion of the drop along the stable flow direction 
towards the stagnation point. Leal's research group at the Chemical Engineering 
Department at UC Santa Barbara had worked since that with the FRM apparatus, mainly in 
the pure elongational flow regimen. In the case of the Parallel Band Apparatus, Birkhoffer 
(2005) proposes computer-controlled flow cell based on the PBA using a digital 
proportional-integral-differential controller. 

We present a nonlinear control applied to study the rheology of a drop in an elongational 
flow field with vorticity. Large deformations on fluid particles such as drops occur typically 
in regions containing saddle points. However, the kinematics of a saddle-point flow does 
not allow for long observation times of the deformed drop due to the outgoing streamlines, 
which advect the drop away from the flow region of interest. Thus, it is necessary to 
accurately control the centroid of the drop to study its desired rheology. A suitable control 
mechanism is essential to maintain the drop position under known flow conditions for times 
that are long compared to the intrinsic time scales of the dynamics. 

In this work, the control mechanism of the position of a drop around the stagnation point of 
the flow generated by a co-rotating Two-Roll Mill is presented. This control is based on the 
Poincare-Bendixson theorem for two-dimensional ordinary differential equations. Namely, 
when a particle moves within a closed region containing a saddle point inside and the 
vector field of the equation points inwards at the boundary of the region, the particle 
undergoes a stable attractive periodic motion. We show that given a prescribed tolerance 
region, around the unstable stagnation point, an incoming flow can always be generated 
when the center of mass of the drop reaches the boundary of the tolerance region. This 
perturbed flow is produced by adjusting the angular velocity of the cylinders, calculated 
using an analytical solution for the flow, without significant change in the flow parameters. 
This gives the time dependent analogue of the Poincare-Bendixson situation just described. 
The drop is controlled in a perturbed attracting periodic trajectory around the saddle point, 
while being confined to a prescribed tolerance area. This mechanism is very different from 
the one used for the proportional control which modifies the unstable nature of the saddle 
point by adjusting the angular velocities of the cylinder to project the motion along the 
stable direction only. In the proposed control the effect of the unstable direction combined 
with the flow readjustment produces the periodic motion. 

A nonlinear control strategy, based on the Poincare-Bendixson theory, is proposed and 
studied. In essence, the control generates the planar motion of a particle (the centroid in this 
case) to ensure a periodic motion of the drop centroid inside a prescribed area around the 
saddle point. In addition, a numerical study and an experimental situation are presented to 
illustrate the effect of the control on the drop motion. The implementation of the nonlinear 
control is within a closed region containing the saddle point, with the velocity field pointing 
inward at every point on the boundary, while the particle undergoes a stable attractive 
periodic motion. Thus, given a prescribed tolerance region around the stagnation point, it is 
always possible to generate a controlled incoming flow whenever the center of mass of the 
liquid drop reaches the boundary of the tolerance region to force the centroid back into the 
prescribed tolerance region. 
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The proposed control scheme is capable of keeping the values of the global flow parameters 
within a small tolerance, and redirecting the liquid-drop centroid toward the stagnation 
point on a time scale that is much shorter than that of the evolution, with minimal impact 
upon the liquid-drop dynamics. 

We also performed detailed studies of the sensitivity of the control to imperfections in the 
shape of the geometry that generates the flow, and the variation of the velocity in the 
servomotors which control the position of the stagnation point. Also a detailed comparison 
between numeric and experimental data is also presented. This provides a complete study of 
the two dimensional situation. Finally, we present some open questions both in the modeling 
and in the theory. In particular in the possibility of extending the current results to drops of 
complex fluids such as viscoelastic drops, vesicles, capsules, and other immersed objects. 

Given that there are no detailed studies available on the influence of the control scheme 
upon the drop's forms, we study numerically the influence of this control on the motion of a 
two dimensional drop by solving the Stokes equations in a container subjected to the 
appropriate boundary conditions on the cylinders and the free surface of the drop. These 
equations are solved with the Boundary Element Method for a variety of flows and drop 
parameters in order to study the perturbation effects introduced by the application of the 
control scheme and to provide the appropriate parameters for future experimental studies. 
In particular, an easy to implement control scheme is an essential tool prior to undertakings 
any experimental studies of the drop's dynamics -under elongational flows with vorticity at 
small Capillary numbers-, for large deformation of drops capsules and other objects, as well 
as for break up and coalescence of embedded objects. 

The proposed method is simple. As the drop evolves under flow conditions, its centroid is 
tracked. When the drop drifts away and its centroid overtakes the prescribed domain about 
the nominal stagnation point, the flow is modified by adjusting the angular velocity of the 
cylinders according to the values obtained from the approximate solution for flows 
generated by TRMs. Essentially, by adjustment the angular velocities of the cylinders, the 
outgoing streamline environment is change into an incoming one, reversing the direction of 
motion of the drop, which is now towards the nominal stagnation point along a stable 
direction. The reversal of direction does not alter significantly the deformation rates applied 
upon the drop; thus, the drop's dynamics is essentially undisturbed. The process is repeated 
as needed and the drop is confined for long times under steady and known conditions. 

This study is complementary to the previous one (Bentley 1986a, Bentley 1986b) because it 
allows detailed studies of the time evolution of the drop's parameters -mainly its 
deformation and the orientation. In contrast to FRMs, the numerical results presented allow 
us to study the dynamics of embedded objects under the nominal flow conditions when the 
drop remains at the stagnation point, as well as the drop in the "controlled flow" with the 
corresponding parameters. The study was carried out both under several viscosity ratios 
and various geometric setups, assessing the robustness of the proposed method and the 
very small influence it has on the drop behavior. The influence of the control scheme on the 
drop's parameters is small with respect to the nominal flow (around 1%). As well, the 
proposed control scheme is capable of relocating the stagnation point on a time scale much 
shorter than the time scale of the drop's evolution. Moreover, this control would remain 
effective during times much longer than the internal time scales for the drop evolution. 



40 Applications of Nonlinear Control 

The fact that a two-dimensional numerical model for the drop is used does not compromise 
the results here presented, mainly because the control scheme should be applicable for small 
deformation of drops or for drops with very small trajectories about the stagnation point. 
For highly elongated drops, capsules and other complex objects, the simulated values of the 
drop's shape may differ from the experimental values, but the applicability of the control 
scheme should be equally robust. 

2. Formulation of the control problem 

As already demonstrated by Bentley (1986a), the only way to maintain fixed the position of 
the drop with respect to the flow field is by changing the location of the stagnation point via 
adjustments of the angular velocities of the cylinders, with the constraint that these changes 
must avoid significant modifications of the flow field. Consequently, a useful control 
scheme for flows by TRMs or FRMs has to maintain the drop as close as possible to the 
stagnation point for a sufficiently long time, making possible studies of the drop dynamics. 
From now on, the selected flow field conditions of a TRM are called nominal, and its 
properties such as the shear rate, flow-type parameter and the position of the stagnation 
point will be denoted by the subscript Norn. 

Figure 1 shows the streamlines around the stagnation point of the unperturbed flow field 
generated by a co-rotating TRM. When a liquid or rigid particle is placed around the 
unstable stagnation point, the particle drifts in the direction of the outgoing streamlines. The 
objective of the control is to maintain a drop about the stagnation point of the nominal flow 
conditions. To construct the control scheme, the trajectory of the centroid of the drop is 
analyzed as the solution of a two dimensional dynamical system. In this case we assume 
that the centroid is away but near the unstable saddle point. Now if the vector field is 
oncoming on the boundary of a box surrounding the unstable stagnation point the system 




Fig. 1, Streamlines generated by a co-rotating Two-Roll Mill with different radii, showing 
the stagnation point in the gap between the rollers, g. The position of the stagnation point 
along the vertical can be moved changing the angular velocities of the rollers. The value of 
the flow-type parameter is a function of the position of the stagnation point. 
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will settle into a periodic orbit inside the box provided the vector field is time independent. 
When there exists time dependence, a bounded motion (which can be periodic or quasi- 
periodic) is produced inside the box. This motion is equally robust as the time independent 
case, confining the centroid of the drop to any prescribed region provided the appropriate 
incoming flow can be produced at the boundary. 

In Fig. 2a, a rectangle is shown about the stagnation point of the flow field. The boundaries 
of this rectangle serve as the limits where the position of the center of mass of the fluid 
particle is allowed to stay at the nominal flow conditions. Fig. 2b shows a detailed sketch of 
the tolerance domain PQRS above and next to the nominal stagnation point. The dark 
streamlines correspond to the nominal flow field with outgoing flow exiting through the 
side QOr and entering at the side ROl of the tolerance area. The dashed flow lines which 
correspond to the shifted stagnation point yss with the flow field essentially reversed: the 
flow enters the PQ side and exists via RS. Also the flow lines which correspond to -yss 
behave in a symmetric manner relative to the symmetric tolerance area. 
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Fig. 2. Streamlines around the stagnation point of the unperturbed flow field generated by a 
TRM. 



When a fluid particle is placed around the stagnation point and the flow is started, the 
particle's centroid drifts away. Assume the centroid is initially at the position A -at time t = 
0- located inside the tolerance area shown in Fig. 3. In this position, the flow corresponds to 
the nominal conditions, and the centroid is subsequently advected along the outgoing 
direction reaching B at t = fon when the control is applied. The control scheme effect is to 
displace the stagnation point to yss, switching the flow field at B to one towards the nominal 
stagnation point. As a result, the centroid follows the flow lines along the path BC, arriving 
at C at time t = f ff- At f ff the flow is reset to the nominal conditions and the stagnation point 
is moved back to the j/jvom position; thus, the centroid follows the path CD. At D the situation 
is repeated but now shifting the stagnation point to - yss until the centroid reaches E where 
the stagnation point is shifted back to the nominal value and the centroid moves towards F 
where the process is repeated. It is to be noted that the transit time spent along paths AB, 
CD, E¥ is much longer that the transit time along sections BC and DE because of the small 
eigenvalues associated with the corresponding incoming directions which produce the 
motion shown. The slope of /,„ can be modified in order to adjust the instant when the flow 
is reset to the nominal conditions. 
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Fig. 3. Trajectory ABCDEF followed by the centroid of a fluid particle in the controlled flow 
field. The control area is shown in grey. The nominal flow corresponds to the darker 
continuous lines, and the dashed lines show the relative displacement of the flow field 
during the controlled portion of the cycle. The angle between the incoming and outgoing 
streamlines at the nominal stagnation point 0Nom and the angle at the corrective flow 6 SS 
have essentially the same values. 

The purpose of the implemented control scheme is to produce always an incoming flow for 
the drop at the boundary of the tolerance area. Thus the center of mass is effectively moved 
as a dynamical system with an unstable rest point (the stagnation point) but with an 
incoming vector field in an area surrounding the box. This arrangement guarantees the 
existence of an incoming vector field that provides a fixed periodic solution by the Poincare- 
Bendixson theory; see Ross (1984). In this case since the incoming vector field is time 
dependent, a bounded trajectory is obtained that is approximately periodic. This nonlinear 
procedure of balancing the repulsion at the critical point with the correction of the boundary 
of the tolerance region always produces a very robust bounded trajectory inside any 
prescribed area. 

All displacements of the stagnation point are assumed to be carried out on a time scale small 
compared to the dynamics of the drop. In the theoretical description above, both the 
centroid of the drop and the streamlines are adjusted instantaneously. For a laboratory 
experiment this will not be the case: the exact position of the centroid is determined after 
processing the flow field images, determining the centroid position and then the driving 
motors modify accordingly the flow field within a finite response time. The latter time lags 
are taken into account in order to calculate the effect of the control in experimental 
situations. 

The relevant times involved are x t associated to the velocity of the video system -fps- and 
the time of capture and processing of all images, the finite response time i 2 of the cylinders 
to readjust their velocity as a consequence of the control, mainly due to the inertia of the 
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mechanical system, and the time of response i 3 of the fluid around the drop to the 
adjustment in the velocity of the TRM. The time t 3 of adjustment of the flow can be 
estimated as the diffusion time t 3 =1 2 /v based on the gap g between the cylinders and the 
kinematic viscosity v. The total response time Tj + i 2 + t 3 = t c must be smaller than the 
characteristic time id of the internal motion of the drop which is a function of the Capillary 
number and the viscosity ratio. In the present work, t x + i> » t 3 given the material 
properties of the fluids and geometry of the setup; Stokes solution implies such time scales 
as well. 

To determine the adjusted velocity field -i.e., reposition the stagnation point position yss (or 
-yss) needed to ensure an incoming flow along PQ - a new flow line (assumed straight) is 
calculated, entering at PQ and ending at yss; see Fig. 3. This requirement gives yss as a 
function of the size of the control area. Parametrizing the flow, relating the position of the 
stagnation point with the angular velocities of the rollers, a x and a 2 , maintaining II2D 
constant, the required values of a c x and a c 2 needed to move the position of the stagnation 
point to yss can be calculated. 

The actual control step is modeled as follows: When the centroid reaches the boundary of 
the tolerance area at B in Fig. 3, the angular velocities are readjusted to the calculated values 
o c t and a c 2 according to the ramp function 

CO 1 — ^Nom if ft — ^c\ \ 

co c i(t) = co Nom ,i -I l + tanh[— - — H (1) 

where t c is the time when the centroid of the particle reaches the point B, and i takes the 
values 1 or 2. The following step at point D is calculated in the same manner. It is remarked 
that during the control steps the flow parameter II2D is kept constant while the change in a is 
less than 0.5% for all cases. 

This control is different from that of Bentley and Leal. In detail, this control does not aim to 
stabilize the center of mass of the drop at the stagnation point as in the proportional control. 
The present control takes advantage of the knowledge of the local flow field and balances 
the unstable motion at the stagnation point with a time dependent incoming flow at the 
boundary, giving an effective dynamical system with a periodic or quasi-periodic orbit for 
the centroid. This is achieved by repositioning the velocity field; thus, producing a bounded 
motion inside the tolerance area as a result of a balance of instability and the modified flow. 
In contrast, the control by Bentley and Leal modifies the flow to project the trajectory of the 
centroid along the stable direction at the nominal stagnation point; essentially, a 
proportional method attempts to trace the centroid of the drop back to the stagnation point. 

3. Experimental setup 

The Two-Roll Mill experimental setup is shown in Fig. 4. The setup consists mainly of the 
Two-Roll Mill flow cell with controlled temperature, the driving system, the optical system 
and the interface system. The flow cell consist of the main body, a set of rollers of different 
radii with machining tolerances for cylinder's diameters and gap of less than 5 urn. The 
parallelism and eccentricity of the rollers axes is limited to less than 5 um , for a top to bottom 
distance of 10 cm. The driving system consists of two servomotors Kollmorgen AKM-11B with 
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Fig. 4. Experimental setup. The flow cell with the motors and the optical system are shown. 

their controllers SERCOS Servostar 300. The optical systems is made up with a Navitar 
microscope -with a telecentric objective with a motorized magnification of 12X-, an adapter 
Navitar 1-61390 with a magnification of 2X, and an IEEE 1394 CCD camera, model XCD- 
X700, by Sony set for a capture rate of 15 and 30 fps. The visual field is about 1.2 mm 
lengthwise covered with 1024 by 768 pixels. The optical resolution (based on the real part of 
the Optical Transfer Function) of the full assembly is about 8 um -or better than 64 line 
pairs/ mm. The main system is mounted on a pneumatically levitated optical workstation by 
Newport Research. 

The typical observed response time for the computer interface, power electronics and 
cylinder's inertia is less than 0.01 seconds, for changes of rotational speeds less than 5% of 
the preset values. The flow parameters and the position of the stagnation point are adjusted 
by varying simultaneously the angular velocities of both cylinders a ± and oo 2 keeping IIzd 
constant. 

Figure 5 shows a schematic block diagram of the experimental setup. A drop is initially 
placed near the stagnation point of the flow cell. The optical system captures images of the 
drop and in the computer, the shape and the center of mass is calculated, and the angular 
velocities of the rollers are calculated in order to maintain the drop in the position desired. 
The calculated angular velocities are sent to the motors controllers and a new image of the 
drop is captured and the cycle is repeated. 

The interface system consists of a workstation HP XW4300 with a PCI SERCOS expansion 
card, which communicates with the motor controllers. The control software is programmed 
in Visual C++, in real-time mode. A Graphical User Interface (GUI) is used to provide access 
to the functionality of the application, see Fig. 6. It incorporates two different aspects that 
can be manipulated separately, one concerning to the video and the other concerning to the 
control of the motors. 

For the video aspect, the GUI window has two displays, one is for the video input, that 
shows the frame that is acquired by the camera; and the other display is for the processed 
image, showing the contour of the drop and its center of mass, along with the tolerance area 
(fixed with the slide bar in the window) and the lines that corresponds to the ingoing and 
outgoing streamlines of the nominal stagnation point. 
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Fig. 5. Experimental block diagram. 
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Fig. 6. Graphical User Interface (GUI) of the Two-Roll Mill Experiment. 

For the control of the motors, the GUI has two parts, one is for the manual control and the 
other is for the automatic control. The manual control window is used for the initial 
positioning of the drop around the nominal stagnation point, inside the tolerance area. 

For the control of the motors the window has slide bars that allow controlling the velocities of 
each cylinder as well as its direction of rotation. In order to place the drops in the center of the 
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image, it is necessary to put it in that location; this is done by manipulating its trajectory with 
the controls above mentioned. Once the drop is in the right place (it means the center of mass 
is inside the tolerance area), the automatic control is activated and the experiment goes on. The 
monitoring section in the window allows us to watch the instantaneous velocities of the 
motors, the coordinates of the center of mass and the size of the tolerance area. 

4. Control scheme implementation 

The automatic control for the experiment requires a real-time process. This program consists 
basically of three parts: (i) The image acquisition (it) Image analysis (tit) Calculus and 
adjustment of the velocities of the motors. 

4.1 Image acquisition and analysis 

The images are provided by the CCD camera using the Instrumental and Industrial Digital 
Camera Application Programming Interface (IIDCAPI) by Sony in the C++ program to 
handle the frames provided by the camera. This section consists in two parts, the first one 
take the frame and stored it in a file, the second one creates a list constantly updated in 
order to always have the last image acquired by the camera available for the analysis section 
and in case of some delay in acquiring, have a reservoir of frames. 

To carry out the analysis of the last frame taken, in order to find the center of mass of the 
drop (in pixels) we use the Open Source Computer Vision Library (OpenCV). The principal 
parameter is the threshold value which can be adjusted ranging from to 255 in a gray- 
scale, used to make a binary image. In the binary image, the contour of the drop is 
computed using the Canny algorithm. Then the center of mass is found using the 
corresponding discretized integral. 



4.2 Calculus and adjustment of the velocities of the motors 

Is in this part where the control takes place. Once the program has the coordinates of the 
center of mass of the drop, decides whether is inside the tolerance area. If is not, the 
program calculates and modifies the velocities of the motors depending on the position of 
the center of mass of the drop. Also, in this part, the data of the number of the frame 
processed, the coordinates of the center of mass and the size of the tolerance area are stored 
in a file. The diagram of the control scheme is shown in Fig. 7. 
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Fig. 7. Diagram of the control scheme. 
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5. Experimental results 

5.1 Parameter used in studies of drop deformation 

The parameters that govern the drop deformation and breakup are the ratio of the drop 
viscosity to that of the suspending fluid A M , the tensorial character of Vu, the history of the 
flow and the initial drop shape. 

The Capillary number Ca represents the ratio of flow forces to surface tension, it is given by 

a/ft a|H 2Z ,| 1 /2 J u 1 



Ca ■ 



(2) 



Where II 2D is the second invariant of 2D = Vu + Vu T . The tensorial character of Vu gives the 
flow- type parameter a. For a given elongational flow with vorticity, a values close to 1 
imply an elongation dominated flow, while values close to zero imply a flow with strong 
vorticity; that is a is a measure of the of the strength of the flow causing drops to deform, 
while the vorticity present in the flow induces s rotation of drops and can inhibit drop 
breakup. From the definition of the flow-type parameter a, (see Astarita, 1979), 



1 + a 



IDI 



(3) 



Thus a is given by 



|D||-||W|| 



\\D\\ + 



(4) 



Where W is the objective vorticity tensor which measures the rate of rotation of a material 
point with respect to the rate of deformation's principal axes at that point. 

A dimensionless measure of the magnitude of the drop deformation is needed. This 
parameter is the Taylor Deformation Parameter Dt, defined in terms of the longest and 
shortest semi-axes of the ellipsoidal drop cross section, see Fig 8. 




Fig. 8. Deformation Parameter DT and orientation angle. 



'max 'mm 



+ r„ 



'max ' 'mm 



(5) 
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The orientation angle of the drop is the angle between the longest axis of the drop and the x- 
axis. 

The flow-type parameter value is related with the angles between the incoming and the 
outgoing axes shown in Fig. 3 as 8 Nom — 2 arctan^ a Nom ) and as 6 ss — 2 arctan{^ja ss \ 
In the experiment, the outgoing flow direction, about the stagnation point, is always 
observed responding to small differences of the refractive index of Fluid 1 due to very small 
differences of temperature between fluid volumes around each cylinder. Based upon 
reversal of the flow field -counter and co-rotating directions of the cylinders-, the angle 
between incoming and outgoing flow can be accurately measured within 0.1 degrees inside 
the visual field -less than 1.2 mm lengthwise. The observed angles are correct within ± 0.05 
degrees of the nominal values here presented. This small angular uncertainty implies 
uncertainties for the nominal flow- type parameter of less than ± 0.5%. Thus, the 
experimental flow obtained is an excellent approximation for the theoretical one predicted. 

The exterior fluid (Fluid 1) is a Polydimethylsiloxane oil DMS 25, n = 485 mPa s, with a 
relative density of 0.971. A drop (Fluid 2) of vegetable canola oil, filtered through a 3 um 
pore size. At 23°C, the viscosity of Fluid 2 is n = 72.6 mPa s, with a relative density of 0.917 is 
used. Both liquids have a well defined Newtonian behavior at the interval of shear rate 
values used. The following figures show the effect on the deformation, orientation and 
trajectories of the centroid of a drop due to variations on the /,„ parameter. This parameter 
permit adjusts the control in order to decrease the drifting effects of the r c time, which is a 
nonlinear function of the viscosity of the surrounding fluid. The drop tested had a diameter 
of 1.0mm, Ca = 0.1031. 



Deformation and Orientation of the drop 




30 t [s] 40 




Fig. 9. Deformation, orientation and trajectory of the centroid of the drop, using the 
parameter h n = 20°. The mean deformation is Dt = 0.1039, STD=0.002 and the mean 
orientation angle is 41.8°. 
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Fig. 10. Deformation, orientation and trajectory of the centroid of the drop, using the 
parameter /,„ = 30°. The mean deformation is Dt = 0.1042, STD=0.0021, the mean orientation 
angle is 41.8°. 
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Fig. 11. Deformation, orientation and trajectory of the centroid of the drop, using the 
parameter /,„ = 40°. The mean deformation is Dt = 0.1045, STD=0.0025, and the mean 
orientation angle is 41.8°. 
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Fig. 12. Deformation, orientation and trajectory of the centroid of the drop, using the 
parameter /,,, = 50° . The mean deformation is Dt = 0.1049, STD=0.0027, the mean orientation 
angle is 41.8°. 



Deformation and Orientation of the drop 




Fig. 13. Deformation, orientation and trajectory of the centroid of the drop, using vertical 
limits. The mean deformation is Dx = 0.1039, STD=0.0020, the mean orientation angle is 41.8°. 
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6. Numerical results and comparisons with the experimental results 

In order to study numerically the evolution of the drop on the proposed control strategy, we 
use the results of Reyes et al. (2011). We just outline the main steps to introduce the 
extensions used in this Chapter to study the effect of noise and the imperfections present in 
the experimental situation. 

We thus proceed as in Reyes et al. (2011) advancing the drop boundary using the equation 

dx 

— = (v, • n)n (6) 

at 

where the velocity v s in Eq. 6 is determined solving a boundary integral equation on the 
drop surface (See Pozrikidis, 1992). 

The control is implemented as follows: From the new drop surface obtained from Eq. 6, we 
determined the center of mass of the drop using the same integral of the experiments. Once 
this is determined, we verified if it falls inside the tolerance region. When it is out of the 
selected boundary, we applied the correcting flow obtained adjusting the angular velocities 
using Eq. 1. 

We considered the effect of noise and imperfections as follows. In the first place, we noted 
from the experiments a systematic variation of the angular velocity due to the geometrical 
imperfections of the cylinder. This was fitted with a single harmonic function with 
frequency and amplitude determined from the experimental values for the observed flow 
without drop. The random noise was taken to be white noise. With these new elements, we 
used the same code described in Reyes et al. (2011) to calculate the trajectories of the drop's 
center of mass, the deformation parameter and the angle of alignment in order to compare 
with the experimental results of the previous Section. 

It is to be noted that the solution of the system are very sensitive in detail to the initial 
conditions. However, the broad features which depend on the limit cycle nature of the 
motion for the center of mass are very robust. 

Because of this reason, we start the numerical solution with initial conditions for the drop's 
center of mass which are taken from the experimental data when the initial rapid transients 
have subsided. Moreover, the numerical flow is started at nominal values since the 
inhomogeneities presented prevent the analytical construction of the initial flow. 

With this, we expect a very good agreement between the numerical and experimental values 
of the deformation and orientation. This is shown in Figures 14-18. Although those figures 
were generated with different initial conditions, the broad behavior is similar to the 
experimental data. The trajectories are also compared. We see good agreement in the broad 
features, in particular, when the control is operational. 

The experimental data shows larger excursions from the nominal stagnation point. These 
are due to the mismatch between the commercial worm gear and worm mechanism used to 
reduce the angular velocity of the motors and transmit the motion to the rollers. 

In Fig. 19 we display the X and Y component of the motion for the center of mass, the blue 
lines shows the experimental values and the red lines the numerical solution. The 
comparison is good considering the mismatch between the initial flows up to t = 10 s. 
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Fig. 14. Trajectory of the centroid of a drop, using vertical limits. The insert graph shows the 
deformation and the orientation angle. The mean deformation is Dt = 0.1039, 
STD = 0.000972 with a mean orientation angle of 41.8°, STD = 0.07626. 
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Fig. 15. Trajectory of the centroid of a drop, using the parameter /,„ same as the incoming 
axis (19.435 °). The insert graph shows the deformation and the orientation angle. The mean 
deformation is D r = 0.10388, STD=0.00104017 with a mean orientation angle of 41.801°, 
STD = 0.0764. 
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Fig. 16. Trajectory of the centroid of a drop, using the parameter /,„ = 30 °. The insert graph 
shows the deformation and the orientation angle. The mean deformation is Dr = 0.103915, 
STD=0.0009935 with a mean orientation angle of 41.7993°, STD = 0.0784447. 




Fig. 17. Trajectory of the centroid of a drop, using the parameter /,„ = 40 °. The insert graph 
shows the deformation and the orientation angle. The mean deformation is Dr = 0.103883, 
STD=0. 00092441 with a mean orientation angle of 41.7992°, STD = 0.0735669. 
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Fig. 18. Trajectory of the centroid of a drop, using the parameter /,„ = 50°. The insert graph 
shows the deformation and the orientation angle. The mean deformation is Dr = 0.103934, 
STD=0. 00096193 with a mean orientation angle of 41.7972°, STD = 0.0733545. 




Y coordinate 




Fig. 19. Experimental and numerical comparisons of the X and Y coordinates of the center of 
mass of a drop subjected to the control. The blue lines are for the experimental data whereas 
the red lines are for the numerical data. 
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Beyond this time, the experimental data shows larger excursions. These are due to the 
mismatch between the commercial worm gear and worm mechanism used to reduce the 
angular velocity of the motors and transmit the motion to the rollers. 

It is to be noted the remarkable agreement in the actual drop parameters which are the ones 
of interest. 

7. Conclusions 

We have shown that it is possible to maintain the position of a drop about the unstable 
stagnation point of the flow field generated by a TRM setup. This scheme is of the upmost 
importance for studies of drops in elongational flow with significant vorticity, and of 
relevance because its space of parameters is not accessible to FRM flows previously studied 
since Bentley (1986a). Indeed, TRMs expand the family of 2D elongational flows amenable 
with a four roll mill, albeit the former flows carry amounts of vorticity similar to that of 
simple shear flows. 

But because the TRMs configuration can only displace the stagnation point along the line 
between the cylinder axes, the control scheme developed for FRMs or PBAs cannot be used 
for studies of drop's dynamics in elongational flows of the TRM type. For TRMs the control 
scheme is based upon features of Poincare-Bendixson limiting cycles. However, these limit 
cycles do not imply that the wandering trajectory of the drop has to be contained in a very 
tight region about the nominal stagnation point. The control scheme appears capable of 
working for a tolerance region with the cases tested experimentally. 

The control implemented in the experimental has shown to be successful. At this point, 
several complications have been seen in the implementation: nevertheless the control 
scheme is robust enough to keep the drop inside a region where the parameters of interest 
have a low variation, for a long times, enough to have reliable measures of the relevant 
parameters. The figures 9-13 shown that despite the trajectory of the drop, the parameters of 
interest in drop dynamics (Dt and orientation angle) are the same. 

It is important to mention that even when the comparisons are made just for only one flow- 
type parameter, the results shown that it is reasonable to expect the same when we will use a 
different flow-type parameter, i.e. a different geometry. 

Important differences have been observed between the experimental and numerical 
trajectories, this disagreement is due to mechanical imperfections, and could be fixed in the 
future by using an system motion transmission with better precision. 
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1. Introduction 

Several diesel-generator sets usually operate on parallel connection in ship power system, 
which has altitudinal nonlinearity. When operating point of system changes, its dynamic 
property will change markedly. The oscillation phenomenon of ship power system which is 
acyclic, random and gusty or paroxysmal will occur on light load working condition, it can 
result in system sectionalizing when it is serious, this phenomenon is called chaos. Chaos is 
a very complicated phenomenon which is generated by the interaction of each parameter in 
the nonlinear system. When it appears in ship power system, following the continuous and 
random oscillation of system operating parameter, which endangers operation security of 
system seriously, it must be prevented and eliminated effectively in the system. In order to 
analyze the chaos phenomenon of ship power system, the nonlinear mathematical model of 
two diesel-generator sets on parallel connection is built in this paper, which reflects the 
variation law of ship power system. Then the light load working condition of two diesel- 
generator sets on parallel connection in ship power station is analyzed by using Lyapunov 
index method on the base of this, seeking the generating mechanism of chaos. A nonlinear 
robust synthetic controller is designed which is based on the nonlinear mathematical model 
of diesel-generator set, then a nonlinear robust synthetic control law is developed for the 
diesel-generator set, it will be applied to control the chaos phenomena, thus providing 
desirable stability for ship power system. 

2. Mathematic model of diesel-generator set on parallel connection 

The mathematical model of diesel-generator set include the mathematical model of 
electromechanical transient process and electromagnetism transient process, first building 
the mathematical model of electromechanical transient process, then the mathematical 



Foundation item: Supported by the National Natural Science Foundation of China under Grant 
No. 60774072; Fundamental Research Funds for the Central Universities of China under Grant 
No.HEUCFT1005. 



58 



Applications of Nonlinear Control 



model of electromagnetism transient process, the mathematical model of one diesel- 
generator set is get on the base of this. 

The mathematical model of diesel-generator set electromechanical transient process 
describes the motion law of diesel-generator set! 1 " 3 !, reflecting the dynamic change process of 
power angle and angular speed, its expression is 



ds 

At ' 


--(0-1, 






da> 


T b 1 C 2 r 

= — — a> + Cj + — — L- 

?>o T «®o T «®o 


1 W . s 

■ — sin o - 

T a°>0 X 'd 


1 U 2 X' d -X q 

7>0 2 X 'd X q 



sin 28. 



(1) 



In the equation: S is power angle of diesel-generator set, a> electric angular speed of diesel- 
generator set, U terminal voltage of generator stator winding, E' q-axis transient electric 
potential, X reactance of each winding, L output axis displacement of diesel engine 
govonor actuator, o =100nrad/s, T a ,T b ,c 1 ,c 2 constants, S,L real values, other variables 
per unit values. 

From Eq.(l) we know, this equation has nonlinear feature. 

The mathematical model of diesel-generator set electromagnetism transient process include 
stator voltage balance equation of generator and electromagnetism transient equation of 
field winding! 3 !, omitting the effect of damping winding, its expression is 



di 
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d j 
1 d' 



U q =-RI q -coX' d I d 
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■pd+tf. 
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(2) 



In the equation: U is terminal voltage of stator winding, U d and U„ d-axis and q-axis 
component of stator winding terminal voltage, R resistance of stator winding, X reactance 
of each winding, I current of each winding, T time constant of each winding, E' q-axis 
transient electric potential, E« voltage of exciting winding. 

Since E' is not easy to measure, we select the terminal voltage of stator winding U as state 
varible, only needing change E' of Eq.(2) into U . According to the relation between varible 
data, the following form is set up 



U = E' q -c 3 S 



(3) 



Substituting Eq.(3) into the first item Eq.(2), we have 
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Substituting Eq.(3) into Eq.(l) and combining with Eq.(4), we have 
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We know from the expression of current I d 

I 
Substituting Eq.(3) into Eq.(6), we have 
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Substituting Eq.(7) into Eq.(5), we get 
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(8) 



Eq.(8) is the nonlinear mathematical model of one diesel-generator set, which reflects the 
relationship of interaction and mutual influence among power angle, speed and voltage, 
describing the variation law of three variables more exactly. 

We select d-q axis of first synchronous generator as reference frame, building the 
mathematical model of two diesel-generator sets on parallel connection. Suppose two diesel- 
generator sets have the same power, type and parameters, d,q component of load current is 
I d ,I , power anger difference of two synchronous generators is S 12 , the mathematic model 
of diesel-generator set on parallel connection is 
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In the equation: i '. = 1, 2 , subscript 1 indicates the first diesel-generator set, subscript 2 
indicates the second diesel-generator set. 

Current coupling relation of two diesel-generator sets is 



<?i 



COS fi- 



ll 



sin <J 



sin S 12 cos 8- 



12 



<1 2 . 



(10) 



Eq.(10) describes the current distribution relation after two diesel-generator sets enter into 
parallel connection. 



Voltage coupling relation of two diesel-generator sets is 



■7 1 



cos S^ 
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(ii) 



Eq.(ll) describes the voltage restriction relation after two diesel-generator sets enter into 
parallel connection. 

Eq.(9), Eq.(10) and Eq.(ll) are the nonlinear mathematical model of two diesel-generator sets 
on parallel connection, which reflects the relationship of interaction and mutual influence 
between two diesel-generator sets, describing the variation law of power angle, speed and 
voltage on two diesel-generator sets exactly. 

3. Chaos oscillation analysis of diesel-generator set on parallel connection 

This paper will make research on the stability of ship power system by the variation law of 
power angle, speed and voltage on diesel-generator set. Due to the power transmission 
between the diesel-generator sets on parallel connection, it is apt to produce power 
oscillation. Power oscillation is the dynamic process of power of diesel-generator set 
regulating repeatedly under the effect of some periodic interference. 

Power oscillation is a chaos oscillation in nature from the point of view of chaos. The 
fundamental feature of chaotic motion is highly sensitive to initial condition, the track which 
is produced by two initial values which is very near each other, will separate according to 
index pattern as time elapses, Lyapunov index is the quantity which describes the 
phenomenon. Distinguishing methods of time series chaotic character include fix quantity 
analysis and ocular analysis, first making numerical analysis for Lyapunov index and 
judging the condition of chaos emerging, then determining if the chaos exists or not under 
this condition by the method of ocular analysis. The methods of ocular analysis include time 
course method, phase path chart method, strobe sampling method, Poincare cross section 
method and power spectrum method. The methods of calculating Lyapunov index include 
definition method, Wolf method and Jacobian method, Jacobian method is a method of 
calculating Lyapunov index which develops in real application. This paper will use Jacobian 
method to calculate Lyapunov index. 

Considering following differential equation system 



x = F(x) 



(12) 
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dx 
In the equation: x = — , x e R m . The evolution of tangent vector e of dot x(t) in the tangent 

At 

space can be expressed by the equation as follow 

e = T(x(t))e,T = ^- (13) 

dx 

In the equation: T is Jacobian matrix of F . The solution of equation can be expressed as 

e(t) = U(t,e(0)) (14) 

In the equation U : e(0) \— > e(t) is mapping of linear operator. The asymptotic behavior of 
mapping U can be described by index as 

A(x(0), e (0)) = hmilnM (15) 

f-*=° t k(0) 

So, the Lyapunov index of system (12) can be formulated as the mean of above repeat 
process 

k-xakAtfi \\e(iAt)\\ 

j-i II - -II ( 16 ) 



-lim 1 ln H e((fc + 1)Af) i i e(fcAf) ll H e(2Af) ll 

k^cokAt \\e(kAt)\\ \\e((k-l)M)\\'" \\e(At)\\ 

For a phase space of n dimension, there will be n Lyapunov index, arranging them 
according to the order from big to small, supposing (A 1 >-l 2 >-">A n ), Xy is called 
maximum Lyapunov index. Generally speaking, having negative Lyapunov index 
corresponds with contracting direction, the tracks which are near are stable in the part, 
corresponding periodic motion. The positive Lyapunov index indicates that the tracks 
which are near separate by index, the strange attractor is formed in phase space, the 
Lyapunov index X is bigger, the chaotic nature of system is stronger, vice versa. For a phase 
space of n dimension, the maximum Lyapunov index is whether bigger than or not is the 
basis of judging the system if has chaos oscillation or not. 

Computing Lyapunov index of two diesel-generator sets on parallel running with light load 
separately, two diesel-generator sets all use conventional controllers. Fig.l and Fig.2 give the 
phase diagram of power angle, speed and voltage of two diesel-generator sets on parallel 
running with 12.5% load separately. 

Two diesel-generator sets run for 100 seconds on parallel connection with 12.5% load, the 
initial value of No.l diesel-generator set is: (S,a),U) = (0.1017,1.0662,0.9762), Lyapunov 
index is: ^=0.076789, /t, = 0.035235 , -i, = -0.197558; the initial value of No.2 diesel- 
generator set is: (S,a>,U) = (0.1022,1.0662,0.9762) , Lyapunov index is: ^=0.076806, 
^2 = 0.035230 , X i = -0.197571 . 

Fig.3 and Fig.4 give the phase diagram of power angle, speed and voltage of two diesel- 
generator sets on parallel running with 25% load separately. 
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Fig. 1. Phase diagram of No.l diesel-generator set when two sets load 12.5% on parallel 
connection 
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Fig. 2. Phase diagram of No.2 diesel-generator set when two sets load 12.5% on parallel 
connection 
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Fig. 3. Phase diagram of No.l diesel-generator set when two sets load 25% on parallel 
connection 
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Fig. 4. Phase diagram of No. 2 diesel-generator set when two sets load 25% on parallel 
connection 



Two diesel-generator sets run for 100 seconds on parallel connection with 25% load, the 
initial value of No.l diesel-generator set is: (S,a>,U) = (0.1812,1.0607,0.9686) , Lyapunov 
index is: ^ = 0.079251 , ^=0.034251, 2, = -0.199955 ;the initial value of No.2 diesel- 
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generator set is: {S,a?,U) =(0.1820,1.0607,0. 

^2 = 0.034953 , A, = -0.199742 . 



Lyapunov index is: /^ = 0.078311 , 



From Fig.l to Fig.4 we can see, all maximum Lyapunov indexes of the system are greater than 
0, showing that the system exists chaotic phenomenon. Two diesel-generator sets with light 
load on parallel connection, enter into chaotic state after running a length of time, their specific 
expression are the oscillation of power angle and speed. Two diesel-generator sets on parallel 
connection load the lighter, the oscillation of power angle and speed is severer. The oscillation 
of power angle means the oscillation of power, the reason is the nonlinearity of ship power 
system and the power transmission between the two diesel-generator sets. The controller in 
this paper is proportional controller, which is a linear controller. It can't control the nonlinear 
character of ship power system, it can't average the load in parallel operation control, there is a 
power angular difference between the diesel-generator sets, thus engendering the power 
transmission between the sets, which results in the happening of chaotic phenomenon. 

Fig.5 and Fig. 6 give the phase diagram of power angle, speed and voltage of two diesel- 
generator sets on parallel running with 25% load plus periodicity load separately. The 
periodicity load usually appears in ship power system, it is widespread. 

Two diesel-generator sets run for 100 seconds on parallel connection with 25% load increasing 



periodicity load O.Olsint, the 
(S,a,U) = (0.1812,1.0607,0.9686) , 
I3 = -0.191497 ; the initial 
(S,a,U) = (0.1820,1.0607,0.9686) , 
I3 = -0.189995 . 



initial value of No.l diesel-generator set is: 

Lyapunov index is: A 1 = 0.073257 , X, = 0.031824 , 

value of No.2 diesel-generator set is: 

Lyapunov index is: A 1 = 0.073393 , X, = 0.030161 , 
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Fig. 5. Phase diagram of No.l diesel-generator set when two sets increase periodicity load 
meanwhile load 25% on parallel connection 
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Fig. 6. Phase diagram of No. 2 diesel-generator set when two sets increase periodicity load 
meanwhile load 25% on parallel connection 

From Fig. 5 and Fig. 6 we can see, all maximum Lyapunov indexes of the system are greater 
than 0, showing that the system exists chaotic phenomenon. The oscillation of power angle 
and speed in two diesel-generator sets is severer than the state of not increasing periodicity 
load. Because the periodicity load is nonlinear, increasing periodicity load intensifies the 
nonlinearity of system, thus aggravating the power oscillation of systemi 4 - 6 I. 

The computer simulation results show that it exists chaotic oscillation phenomenon when 
two diesel-generator sets run on parallel connection with light load. Primary cause of 
emerging this phenomenon is the nonlinearity of ship power system, minor cause is the 
power transmission between the two diesel-generator sets. Besides this, using conventional 
linear controller is also a key factor of generating the chaotic oscillation of system. Only 
using nonlinear controller, making the nonlinear characteristic of ship power system offset 
and compensate, can we solve the problem of system chaotic oscillation fundamentally. The 
chaotic oscillation phenomenon is transition state between stable state and unstable state, it 
must be prevented in order to ensure the stability of system. 



4. Design of nonlinear robust synthetic controller 

Mixed H-two/H-infinity control theory is a robust control theory that has speed 
development from the eighties of 20 century, which can solve the problem of robust stability 
and robust performance! 7 - 10 ]. Because diesel-generator set control system is a nonlinear 
control system, using the method of direct feedback linearization to linearize the nonlinear 
system, the state feedback controller is designed for linearization system using mixed H- 
two/H-infinity control theory, thus acquiring nonlinear robust control law in order to reach 
the purpose of restraining the chaotic phenomenon of ship power system, improving the 
stability of ship power system. 
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Because of coupling action between speed and voltage, the nonlinear robust synthetic 
controller is designed for diesel-generator set in order to control speed and voltage 
synthetically, making the both interaction in minimum range, thus improving the stability of 
frequency and voltage in ship power system. 

The principle diagram of diesel-generator set synthetic control system based on nonlinear 
robust synthetic controller is shown in Fig. 7. The diesel-generator set synthetic control 
system is made up of diesel engine, generator, nonlinear robust synthetic controller, 
actuator, oil feeding mechanism and exciter. Nonlinear robust synthetic controller is made 
up of two parts, one is nonlinear H-two/H-infinity speed controller, another is nonlinear H- 
two/H-infinity voltage controller. 



The differential equation of actuator is 



dL 

"dl 



r, r. 



(17) 



The differential equation of exciter is 



dt 



E R+ K 1 
T, T, 



(18) 
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Fig. 7. Principle diagram of diesel-generator set synthetic control system 

First step, we design nonlinear H-two/H-infinity speed controller. 

Combining Eq.(l) with Eq.(17), we can get the nonliear mathemaical model of diesel engine 
speed regulation system. 
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(19) 



Since Eq.(19) has nonlinear feature, the method of direct feedback linearization is used to 
linearize Eq.(19). Order x x = S , x 2 = a> - 1 , 
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In the equation: ijW is the undesired signal which is assumed for using H-two/H-infinity 
control method, including equivalence disturbance which is generated by disturbance 
torque and modeling error. 

Assigning virtual controlled variable 
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So Eq.(20) can be changed as 



x 2 = x 3 + A^w 
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Eq.(22) can be written as matrix form 



x = Ax + BjW + B 2 v 



(23) 



In the equation: 
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Defining the evaluation signal of dynamic performance as 
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D, 



,C 1 ,C 2 ,D n ,D 12 /D2i /D22 are weighting matrix, </,■■ >0(z'=l,2; j=l,2,3)and r, >0 



J 



(i=l,2) weighting coefficient. We can select optimal performance combination through 
changing weighting coefficient, including stability of ship power system, frequency 
regulation precision and low energy loss of speed regulation system. 

For the control system made up of Eq.(23) and Eq.(24), requring design a controller F, 
making the closed loop system asymptotically stable, moreover H-infinity norm of closed 
loop transfer function T<xj(s) from w to z m not more than a given upper bound, in order to 
ensure the closed loop system have robust stability to uncertainty enter from w; meanwhile 
making H-2 norm of closed loop transfer function T2(s) from w to z 2 as small as possible, so 
as to assure the system performance measured by H-2 norm in a good level, this control 
problem is called H-two/ H-infinity control problem. 

From Eq.(23) and Eq.(24), we can get the augmentation controlled object based on mixed H- 
two/H-infinity control theory 



C 9 



D n D 12 
D 21 D 22 



(25) 



controller F can be solved by corresponding augmentation controlled object P. 
For controlled object P, existing H-two/H-infinity state feedback controller: 
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(26) 



In the equation: F is feedback coefficient, which can be got by using /^-Analysis and 
Synthesis Toolbox in the MATLAB. 

Combining Eq.(21) with Eq.(26), we get 
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That is nonlinear H-two/H-infinity speed control law of diesel-generator set. Substituting 
x 1 ,x 2 ,x 3 into Eq.(27), we can get practical form of nonlinear H-two/H-infinity speed 
control law: 






c 2 Ki 2 X' d X q 



,1 ~ X q 



f 3 sin2S- 



E'UT, 



c 2 K x X' d 



— cos<5(«»-l) + 



U%(X' d -X q ) 

c 2 Klx' d x. 



(28) 



cos2S(m-l) 



Robust Control Research of Chaos 

Phenomenon for Diesel-Generator Set on Parallel Connection 



69 



Second step, we design nonlinear H-two/H-infinity voltage controller. 

Combining Eq.(l), Eq.(4) with Eq.(18), we can get the nonliear mathemaical model of 
synchronous generator voltage regulation system. 
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We select the error of voltage All as state varible, The relation between 17 and All is 



All = 17-17, 



(30) 



In the equation: L7 is initial value of stator winding terminal voltage, its value is 1. 
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substituting Eq.(30) into Eq.(29), we get 



d£ 



f* 



-fd , K 2 



dt 

dAL7_J_ 

a>-\ 



T 2 T 2 



dt 

dJ 

dt ' 
dm _ T b 



E a —_ 



AL7 —S-c 3 6) + (c 3 )- 



X A -X' A 



T. 



1 c, 

-c, + - 



1 



dt T a o> T a o> T a a> T a o> 
Eq.(31) can be written as matrix form 
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Defining the evaluation signal of dynamic performance as 



z' m = C[x' + D' u w' + D' 12 u 2 
z' 2 = C' 2 x' + D' 21 w' + D' 22 u 2 



(33) 



In the equation: C\ 
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C[ , C' 2 , D' n , D' 12 , D' 21 , D' 22 are weighting matrix, q i i >0(j=1,2; j=4r- 



"0 0" 







7)and r ; >0(i=3-10)weighting coefficient. We can select optimal performance combination 
through changing weighting coefficient, including stability of ship power system, voltage 
regulation precision and low energy loss of excitation system. 

From Eq.(32) and Eq.(33), we can get the augmentation controlled object based on mixed H- 
two/H-infinity control theory 



A' 



c; 
c; 



D u D 12 
D 21 D 22 



(34) 



For controlled object P' , existing H-two/H-infinity state feedback controller: 



u 2 =F'x' = [f{ f 2 fi ft]- 



AU 

S 



■-flEfd + fi&U + fiS + flco 



(35) 



That is nonlinear H-two/H-infinity voltage control law of diesel-generator set. 

Third step, we design nonlinear robust synthetic controller. 

Combining Eq.(28) with Eq.(35), we can get nonlinear robust synthetic controller of diesel- 
generator set! 11 - 18 ! 
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(36) 



Eq.(36) considers the coupling function of speed and voltage, which controls both 
synthetically. It can increase dynamic precision of speed and voltage, improving the stability 
of ship power system. 



5. Results of computer simulation 

The key parameters of diesel-generator set control system in this paper are as follow: 

The power of diesel-generator set is 1250kW; the rated speed n =1500r/min; the rotary 
inertia of set /= 71. 822kg m 2 ; the damping coefficient of set D =5.54;the magnetic pole pair 
number of generator p=2; the rated torque of diesel engine 11.9kN m; the maximum troke of 
output axis 10mm. 

The rated voltage of synchronous generator is 390V; the rated current 2310A; the power factor 
0.8; the rated frequency 50Hz; the exciting voltage of exciter 83V; the exciting current 7.7A. 

Designing nonlinear synthetic controller based on mixed H-two/H-infinity control theory, 
assuming disturbance signal coefficient of Eq.(23) d r =0.1, assuming weighting coefficient of 
Eq.(24) q n =0.002, q 12 =0.4, q 13 =0.5, r y =0.01, q 21 =0.002, q 22 =0.4, q 23 =0.5, r 2 =0.01. 



Corresponding matrix are 
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Using LMI toolbox, we get state feedback controller: 

F= [-0.3582 -35.1042 -157.7578] 



(37) 



Assuming weighting coefficient of Eq. (33) q u =0.1, g 15 =2800, q 16 =0.1, q 17 =0.1, q 2i =0.1, 



q 25 =2800, q 26 =0.1, q 27 =0.1, r 3 = r i - 
Corresponding matrix are 
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c[=c 2 
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Using LMI toolbox, we get state feedback controller: 

F = [-0.027 -697.798 -0.023 



-0.025] 



(38) 



Substituting Eq.(37) and Eq.(38) into Eq.(36), we get nonlinear robust synthetic controller of 
diesel-generator set. 

Optimization of weighting function is difficult point of H-two/H-infinity control, needing 
select repeatedly. After each selection, using LMI toolbox to get state feedback coefficient, 
substituting simulation model, making charactrsitc test in order to get best combination 
property index. Fig.8 give the block diagram of the system simulation. 









No.l diesel-generator set 




speed 






















nonlinear 

robust 
synthetic 
controller 
















actuator 




diesel 
engine 


generator 




voltage 














► 
















A 


V 












exciter 






































— ► 


load 














No. 2 diesel-generator set 












speed 


























nonlinear 

robust 
synthetic 
controller 


















actuator 




diesel 

engine 


generator 




voltage 














►! 
















A 


L 












exciter 























































Fig. 8. Block diagram of the system simulation 

In order to test and verify the effect of nonlinear robust synthetic controller for restraining 
the chaotic oscillation of ship power system, Fig.9 and Fig.10 give the phase diagram of 
power angle, speed and voltage of two diesel-generator sets on parallel running with 12.5% 
load separately after using nonlinear robust synthetic controller. The initial values of two 
sets are same as that of generating chaotic oscillation, the running time is also 100 seconds. 
Calculating the Lyapunov index of two diesel-generator sets on parallel operation, we get 
the following results. Lyapunov index of No.l diesel-generator set is: ^=-0.003435, 
/I2 = -0.104389 , /I3 = -0.230524 ; Lyapunov index of No. 2 diesel-generator set is: 
Xy = -0.003442 , ^ = -0.104382 , 4 = -0.230524 . 
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Fig.ll and Fig.12 give the phase diagram of power angle, speed and voltage of two diesel- 
generator sets on parallel running with 25% load separately after using nonlinear robust 
synthetic controller. The initial values of two sets are same as that of generating chaotic 
oscillation, the running time is also 100 seconds. Calculating the Lyapunov index of two 
diesel-generator sets on parallel operation, we get the following results. Lyapunov index of 
No.l diesel-generator set is: A x = -0.003983 , A 2 = -0.103822 , /L, = -0.230553; Lyapunov index 
of No.2 diesel-generator set is: \ = -0.003993 , ^ = -0.103811, A 3 = -0.230554 . 

Fig.13 and Fig.14 give the phase diagram of power angle, speed and voltage of two diesel- 
generator sets on parallel running with 25% load plus periodicity load separately after using 
nonlinear robust synthetic controller. The initial values of two sets are same as that of 
generating chaotic oscillation, the running time is also 100 seconds. Calculating the 
Lyapunov index of two diesel-generator sets on parallel operation, we get the following 
results. Lyapunov index of No.l diesel-generator set is: Ay = -0.004331 , A 2 = -0.103293 , 



i 3 = -0.230741 ; 
A 2 = -0.103283 , 



Lyapunov index 

1 3 = -0.230742 . 



of No.2 diesel-generator set 



A, 



-0.004341 , 



From Fig. 9 to Fig.12 we can see, all maximum Lyapunov indexes of the system are less than 
0, showing that the system doesn't exists chaotic phenomenon and works in a stable range. 
Two diesel-generator sets with light load on parallel connection, run after using nonlinear 
robust synthetic controller, their chaotic phenomenon disappears, power angle, speed and 
voltage run nearby the desired values. It shows that nonlinear robust synthetic controller 
can control the nonlinear character of ship power system effectively and make the nonlinear 
characteristic of ship power system offset and compensate, it can solve the problem of 
system chaotic oscillation fundamentally. 
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Fig. 9. Phase diagram of No.l diesel-generator set when two sets load 12.5% on parallel 
connection after using nonlinear robust synthetic controller 
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Fig. 10. Phase diagram of No. 2 diesel-generator set when two sets load 12.5% on parallel 
connection after using nonlinear robust synthetic controller 
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Fig. 11. Phase diagram of No.l diesel-generator set when two sets load 25% on parallel 
connection after using nonlinear robust synthetic controller 
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Fig. 12. Phase diagram of No.2 diesel-generator set when two sets load 25% on parallel 
connection after using nonlinear robust synthetic controller 

From Fig.13 to Fig.14 we can see, all maximum Lyapunov indexes of the system are also less 
than 0, showing that the system also doesn't exists chaotic phenomenon and works in a 
stable range. Although after adding periodicity load intensifying the nonlinearity of system, 
nonlinear robust synthetic controller can restrain the chaotic phenomenon of ship power 
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Fig. 13. Phase diagram of No.l diesel-generator set when two sets load 25% plus periodicity 
load on parallel connection after using nonlinear robust synthetic controller 
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Fig. 14. Phase diagram of No. 2 diesel-generator set when two sets load 25% plus periodicity 
load on parallel connection after using nonlinear robust synthetic controller 

system, making the nonlinear characteristic of ship power system offset and compensate, 
thus improving the stability of ship power system. 

Because ship power system made up of diesel-generator sets is a nonlinear system, using 
nonlinear robust synthetic controller can offset and compensate the nonlinear characteristic 
of ship power system, it can average the load in parallel operation control, there is no power 
angular difference between the diesel-generator sets, thus avoiding the power transmission 
between the sets, solving the problem of system chaotic oscillation fundamentally, 
improving the voltage and frequency stability of ship power system. The research on 
nonlinear robust synthetic controller of diesel-generator set provides a new control method 
for ship power system, having important practical significance and extensive using 
prospect. 



6. Conclusion 

In order to analyze the chaos phenomenon of ship power system, the nonlinear 
mathematical model of two diesel-generator sets on parallel connection is built in this paper, 
which reflects the relationship of interaction and mutual influence between two sets. The 
light load working condition of two diesel-generator sets on parallel connection in ship 
power station is analyzed by using Lyapunov index method, which proves the presence of 
chaos phenomenon. A nonlinear robust synthetic controller is designed which is based on 
the nonlinear mathematical model of diesel-generator set. In combining the direct feedback 
linearization with robust control theory to design the synthetic controller for the diesel- 
generator set, then a nonlinear robust synthetic control law is developed for the diesel- 
generator set. The computer simulation results show that the nonlinear robust synthetic 
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controller effectively suppresses the chaos phenomenon of ship power system, thus 
providing desirable stability for ship power system. 
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1. Introduction 

The design of robust model reference adaptive control (MRAC) schemes for plants in 
controllable form, comprising unknown varying but bounded coefficients and varying control 
gain has attracted a great deal of research. Many nonlinear systems may be described by the 
controllable form; for instance, second order plants (see (Hong & Yao, 2007), (Hsu et al., 2006), 
(Yao & Tomizuka, 1994), (Jiang & Hill, 1999)) and systems whose nonlinear behavior or part 
of it, is represented by some function approximation technique (cf. Nakanishi et al. (2005), 
(Chen et al., 2008), (Tong et al., 2000), (Huang & Kuo, 2001), (Yousef & Wahba, 2009), (Hsu 
et al., 2006), (Koo, 2001), (Labiod & Guerra, 2007)). 

The state adaptive backstepping (SAB) of (Kanellakopoulos et al., 1991) is a common 
framework for the design of adaptive controllers for plants in controllable form. As is well 
known, a major difficulty in introducing robustness techniques to SAB based schemes is that 
the states z,- and the stabilizing functions must be differentiable to certain extent (see (Yao & 
Tomizuka, 1997), (Yao, 1997), (Ge & Wang, 2003)). 

The robust SAB scheme of (Zhou et al, 2004), (Su et al., 2009), (Feng, Hong, Chen & Su, 
2008) has the advantage that the knowledge on the upper or lower bounds of the plant 
coefficients can be relaxed if the controller is properly designed and the control gain is constant 
or known. The approach is based on the truncation method of (Slotine & Li, 1991), pp. 309. 
The stabilizing functions are smoothed at each i-th step in order to render it differentiable 
enough. The following benefits are obtained: i) the scheme is robust with respect to unknown 
varying but bounded coefficients, ii) upper or lower bounds of the plant coefficients are not 
required to be known, and iii) the tracking error converges to a residual set whose size is 
user-defined. 

The specific case of unknown varying control gain is an important issue, more difficult to 
handle than other unknown varying coefficients. The varying control gain is usually handled 
by means of robustness techniques (cf. (Wang et al., 2004), (Huang & Kuo, 2001), (Bechlioulis 
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& Rovithakis, 2009), (Li et al, 2004)) or the Nussbaum gain technique (cf. (Su et al., 2009), 
(Feng, Hong, Chen & Su, 2008), (Feng et al., 2006), (Ge & Wang, 2003)). The above methods 
are applicable to plants in parametric-pure feedback or controllable form, and with controllers 
that use the SAB or the MRAC as the control framework. 

In (Wang et al., 2004), a system with dead zone in the actuator is considered, assuming that 
both dead zone slopes have the same value. The input term is rewritten as the sum of an input 
term with constant control gain plus a bounded disturbance-like term. The disturbance term 
is rejected by means of a robust technique, based on (Slotine & Li, 1991) pp. 309. Nevertheless, 
this strategy is not valid for different values of the slopes. Other robustness techniques 
comprise a control law with compensating terms and either a projection modification of the 
update law, as in cf. (Huang & Kuo, 2001), or a a modification as in (Bechlioulis & Rovithakis, 
2009), (Li et al., 2004). Nevertheless, some lower or upper bounds of the plant coefficients are 
required to be known. 

The Nussbaum gain technique can relax this requirement, as can be noticed from (Su et al., 
2009), (Feng, Hong, Chen & Su, 2008), (Feng, Su & Hong, 2008). The main drawbacks of 
the Nussbaum gain method are (see (Su et al., 2009), (Feng, Hong, Chen & Su, 2008), (Feng 
et al., 2006), (Ge & Wang, 2003), (Feng et al., 2007), (Feng, Su & Hong, 2008), (Ren et al, 
2008), (Zhang & Ge, 2009), (Du et al., 2010)): i) the upper bound of the transient behavior of 
the tracking error is significantly modified in comparison with that of the disturbance-free 
case: the value of this bound depends on the time integral of terms that comprise Nussbaum 
terms, and ii) the controller involves an additional state, which is necessary to compute the 
Nussbaum function. 

Other drawbacks are: i) the control gain is assumed to be the product of a unknown constant 
and a known function, as in (Tong et al., 2010), (Liu & Tong, 2010), ii) the control gain is 
assumed upper bounded by some unknown constant, as in (Zhang & Ge, 2009), (Du et al., 
2010), (Su et al., 2009), (Feng, Hong, Chen & Su, 2008), (Feng et al., 2006), (Ge & Wang, 2003), 
(Feng et al., 2007), (Feng, Su & Hong, 2008), (Ren et al., 2008), iii) the control gain is assumed 
upper bounded by a known function, as in (Ge & Tee, 2007), (Psillakis, 2010), iv) upper 
or lower bounds of the plant coefficients are required to be known to achieve asymptotic 
convergence of the tracking error to a residual of user-defined size, as in (Ge & Tee, 2007), 
(Chen et al., 2009), (Feng et al., 2006), (Ge & Wang, 2003), (Feng et al., 2007), (Ren et al, 2008), 
(Ge & Tee, 2007), (Tong et al., 2010), (Liu & Tong, 2010), iv) the control or update law involves 
signum type signals, as in (Zhang & Ge, 2009), (Du et al., 2010), (Psillakis, 2010), (Su et al, 
2009), (Feng, Hong, Chen & Su, 2008), (Feng, Su & Hong, 2008). 

Recent adaptive control schemes based on the direct Lyapunov method achieve improved 
transient performance. For instance, L\ adaptive control, with the drawback that the control 
gain is assumed constant, as in (Cao & Hovakimyan, 2006), (Cao & Hovakimyan, 2008a), (Cao 
& Hovakimyan, 2008b), (Dobrokhodov et al., 2008), (Li & Hovakimyan, 2008). 

Other works have the following drawbacks: 

i) The control gain is assumed constant, as in (Zhou et al., 2009), (Wen et al., 2009), (Bashash & 
Jalili, 2009). 

ii) The control gain is assumed upper bounded by some unknown constant, as in (Chen, 2009), 
(Ho et al., 2009) and (Park et al., 2009). 
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iii) The control gain is assumed upper bounded by some known function as in (Bechlioulis & 
Rovithakis, 2009). 

iv) Upper or lower bounds of plant parameters are required to be known to achieve 
asymptotic convergence of the tracking error to a residual set of user-defined size, as in 
(Bashash & Jalili, 2009), (Chen, 2009), (Ho et al, 2009), (Park et al., 2009) and (Bechlioulis 
& Rovithakis, 2009). 

In this chapter, we develop a controller that overcomes the above drawbacks, so that: 

Bi) The upper bound of tracking error transient value does not depend on time integral terms. 

Bii) Additional states are not used in the controller. 

Biii) The control gain is not required to be upper bounded by a constant. 

Biv) The control gain is not required to be bounded by a known function. 

Bv) Upper or lower bounds of the plant parameters are not required to be known. 

Bvi) The control and update laws do not involve signum type signals. 

Bvii) The tracking error converges to a residual set whose size is user-defined. 

We consider systems described by the controllable form model with arbitrary relative degree, 
unknown varying but bounded coefficients and varying control gain. We use the SAB of 
(Kanellakopoulos et al., 1991) as a basic framework for the control design, preserving a simple 
definition of the states resulting from the backstepping procedure. We use the Lyapunov-like 
function method to handle the unknown time varying behavior of the plant parameters. All 
closed loop signals remain bounded so that parameter drifting is prevented. 

The key elements to handle the varying behavior of the control gain are: i) introduce the 
control gain in the term involving the adjusted parameter vector, by means of the inequality 
that relates the control gain and its lower bound, and ii) apply the Young's inequality. 

In current works that deal with plants in controllable form and time varying parameters and 
use the state transformation based on the backstepping procedure, they modify the defined 
states at each step of the state transformation in order to tackle the unknown time varying 
behavior of the plant parameters. Instead of altering the state transformation, we formulate a 
Lyapunov-like function, such that its magnitude and time derivative vanish when the states 
resulting from the state transformation reach a target region. 

The control design and proof of boundedness and convergence properties are simpler in 
comparison to current works that use the Nussbaum gain method. The controller is also 
simpler as it does not introduce additional states that would be necessary to handle the 
unknown time varying control gain. 

The chapter is organized as follows. In section 2 we detail the plant model. In section 3 we 
present the goal of the control design. In section 4 we carry out a state transformation, based 
on the state backstepping procedure. In section 5 we derive the control and update laws. In 
section 6 we prove the boundedness of the closed loop signals. In section 7 we prove the 
convergence of the tracking error e, finally, in section 8 we present an example. 
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2. Problem statement 

In this section we detail the plant and the reference model. Consider the following plant in 
controllable form: 

y(») = 7 J a + bu + d (1) 

where y(t) £ R is the system output, u(t) £ R the input, a a vector of varying entries, J n 
a known vector, b the control gain, and d a disturbance-like term. We make the following 
assumptions: 

Ai) The vector a involves unknown, time varying, bounded entries a.\, • • • , flj, which satisfy: 
\a\\ < \i\, ■ ■ ■ , |fl,| < flj, where }i\, ■ ■ ■ , flj are unknown, positive constants. 

Aii) The entries of the vector 7„ are known linear or nonlinear functions of y, ■ ■ ■ , y'"~ ' . 

Aiii) The terms y ,y, • • • , y are available for measurement. 

Aiv) The term d represents either external disturbances or unknown model terms that satisfy: 

\d\ < Vdfi, (2) 

where }i& is an unknown positive constant, and fy is a known function that depends on y, • • • , 
y( n-1 J . In the case that d is bounded, we have fj = 1. The term d may come from the product 
of a known function gj with an unknown varying but bounded coefficient c„: d = c^g^, 
\Cg\ < \i&, so that /,j = \g&\, where }i& is an unknown positive constant whereas g$ is a known 
function. 

Av) The control gain b satisfies: 

\b\ > b m > 0, b ^ Vf > t (3) 

where b m is an unknown lower bound, and the value of the signum of b is constant and 
known. 

Remark. We recall that fi^, b m , fi\, ■ ■ ■ , flj are unknown constants. In contrast, the values of y, 
■ ■ ■ , y(" -1 ', 7„, /^, sgn(fo) are required to be known. Notice in assumption Av that we do not 
require the control gain b to be upper bounded by any constant. That is a major contribution 
with respect to current works that use the Nussbaum gain method, e.g (Su et al., 2009), (Feng, 
Hong, Chen & Su, 2008), (Feng, Su & Hong, 2008), (Feng et al, 2007), (Ge & Wang, 2003). The 
requirement about the value of the signum of b is a common and acceptable requirement. 

3. Control goal 

Let 



e(t) = y(0-w(0 ( 4 ) 

M i (" — 1) i /r\ 

V d + <hn,n-iy d ' + ••• + a mfi y d = a mfi r (5) 

n e = {e : \e\ < C be } (6) 
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where eit) is the tracking error, y d (t) is the desired output, Cl e is a residual set, r is the reference 

signal. Moreover, « mn _i, ■ ■ ■ , a nlj0 are constant coefficients defined by the user, such that 

the polynomial K(p) is Hurwitz, being K(p) defined as K(p) = p( n > + &m n —\V + ' ' ' + 

«jh /0 . The reference signal r(t) is bounded and user-defined. The constant Q, e is positive and 

user-defined. 

The objective of the MRAC design is to formulate a controller, provided by the plant model 

(1) subject to assumptions Ai to Av, such that: 

i) The tracking error e converges asymptotically to the residual set Cl e . 

ii) The control signals are bounded and do not involve discontinuous signals. 

4. State transformation based on the state backstepping 

In this section we carry out a state transformation by following the steps 0, ■ ■ ■ , n of the 
backstepping procedure. The plant model (1) can be rewritten as follows: 

if = X{ + i, 1 < i < n — 1 (7) 

x„ = fl T 7n(xi, • • • ,x„) + bu + d (8) 

x\ = y, x 2 = y, ■ ■ ■ , x n = y( n ~ > 

The model (7, 8) can be obtained by making 7i = • • • = 7 n -l = in the parametric - pure 
feedback form of (Kanellakopoulos et al., 1991). We use the SAB of (Kanellakopoulos et al., 
1991) as the basic framework for the formulation of the control and update laws. 

We develop the SAB for the plant model (7, 8), and introduce a new robustness technique. 
Since the order of the plant is n, the procedure comprises the steps 0, ■ ■ ■ , n, to be carried out 
in a sequential manner. 

Step 0. We begin by defining the state Zj as the tracking error: 

zi=e = y-y d = x 1 -y d (9) 

Step i (1 < i < n — 1). At each z'-th step, we obtain the dynamics of the state Z; by deriving it 
with respect to time, and using the definitions of i, + j provided by (7). For the sake of clarity, 
we develop the step 1 and then we state a generalization for (1 < i < n — 1). 

For the case / = 1, we differentiate Zj defined in (9) and use the definition of X\ provided by 
(7) with f = l: 

Z\ = xi - y d = x 2 - y d = x 2 + q>\ (10) 

<pi = n(yd) = -yd 

where q>\ is a known function of y d . Equation (10) can be rewritten as: 

Zl = -ClZ!+Z 2 (11) 

z 2 = x 2 + cjZ! - y d (12) 

where C\ > 2 is a positive constant of the user choice. The dynamic equation of z 2 is obtained 
by differentiating it with respect to time. The same procedure must be followed until the step 
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i = n — 1. To state a generalization, we express z,- as: 

z< = x i+i + tpu (13) 

cpi = (p i {z l ,---z l ,y d ,y d ,-- ■ ,yf), (14) 

where cpi is a known scalar term, that is function of Z\, ■ ■ ■ , z,-, y d , y d , ■ ■ ■ , y, . Equation (13) 
can be rewritten as: 

z, = -azi + Zi+i, (15) 

2;+i = x i+ i + tpi + c,-z, (16) 

where c, > 2 is a positive constant of the user choice. At the step i = n — 1 we obtain the 
dynamic equation for z n _\ and the expression for z„ as indicated by (15, 16). 

Remark. Notice that the definition of the states Z; is similar to that of the disturbance free case, so that 
a simple design is preserved. This is due to the following facts: 

i) Disturbance like terms are absent in the dynamics X\, • • • , x n -\ given by (7), so that they are also 
absent in the dynamics Z\, • • • , i n -\, as can be noticed in (13). 

ii) Dead zone functions of the states z,- are not used. 

Step n. We obtain the dynamics of z n by differentiating it with respect to time and using the 
expression of x n provided by (8): 

z n = bu + fl T 7„ + (p n + d (17) 

(pn = fn(zi,--- ,z„,y d ,--- , y { d n) ) (18) 

In) 
where (p n is a known scalar that is function of Z\, ■ ■ ■ , z n , yd, ' ' '/ y\ ■ Notice that the 

disturbance like term d and the control input u appear explicitly at the dynamics of z„, at 

the step n of the procedure. Thus, we have completed the state transformation, which allows 

us to develop the controller. 

5. Control and update laws 

In this section we develop the control and update laws, taking into account the assumptions 
stated in section 2 and the goals of section 3. The key elements of the procedure are: 

i) Incorporate the assumptions Ai and Aiv, concerning the unknown time varying parameter 
a.\, • • • ,Ui and the disturbance like term d. 

ii) Carry out a linear parameterization. 

iii) Express the parameterization in terms of adjustment error and adjusted parameter vector. 

iv) Introduce the control gain b within the adjusted parameter vector. 

v) Formulate the control law. 

vi) Formulate a Lyapunov-like function and find its time derivative. 
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vii) Formulate the update law. 

We begin by rewriting (17) as follows: 

z„ = -c n z n + bu + a T 7„ + cp n + c„z„ + d, (19) 

where c„ > 2 is a positive constant of the user choice. Multiplying (19) by z n , we obtain: 

z„z„ = -c n z\ + bz„u + z„a T j„ +z„((p„ + c„z„) + z„d (20) 

The term z n a j n + z n (q> n + c n z n ) + z n d can be rewritten as follows: 

z n a T 7„ + z„(q>„ + c„z„) + z„d = z„(a T y n + d+ cp n + c n z n ) 

= z «(«[l]?n[l] H ha [/']Tn[j] +d + fn + C n Z„) (21) 

< l z "l (l fl [l]%.[l]l + ■ ■ ■ + \ a [i]7n[j}\ + \ d \ + \<Pn + C„Znl) 

using assumptions Ai and Aiv of section 2 and parameterizing, we obtain: 

z«fl T ?« +z„(ip„ +c n z n ) + z n d 

< l z «l (N7„[i]\ + ■ ■ ■ + P)hu[j]\ + Hdfd + \fn + c„z n \^ (22) 

= VK^\z n \<p T 8 

where 

f = [l7„[i]|/ •" > |7 n [/]l/ Si, Wn + c„z„|] (23) 

e = {i/Vhm) [pi, ■■■ , flj'Fd, i] (24) 

Notice that the entries of the vector 9 are unknown, positive, constant, because bounds of 
the time varying parameters fl, and d have been introduced, according to the properties in 
assumptions Ai and Aiv of section 2. Now, we express (22) in terms of adjustment error and 
adjusted parameter vector: 

Z„a T J n + Z n (cp n + C n Z n ) + Z n d < -\Jb mn \z n \^ Q + \fVmn\Zn\if 9 (25) 

where 



= 6--i==W ■■■,H P m, if (26) 

sjbmn L J 

being 9 an adjusted parameter vector and 9 an adjustment error. Using the property (3) in the 
term \Jb mn \z n \ q> T 9 of (25), we obtain: 

3Clwz : " \b\\z„\\f 0\ (27) 



■be 



where C bvz = (l/2)C 2 be (28) 



86 Applications of Nonlinear Control 

using the Young's inequality (cf. (Royden, 1988) pp. 123), we obtain: 

VKn\z„\\<p T e\ < \c hvz + 3^H^«(<p T e) 2 (29) 

Substituting (29) into (25), we obtain: 

z«n T 7n +z„((p„ + c„z„) + z n d 



< -Vbmn\z„\cp T 6+ \C bm + or— \b\zl(<p 



(30) 



Remark. We have proposed a new method to handle the unknown varying behavior of the control gain 
b, alternative to the current Nussbaum gain method. We parameterized the model term z„a T 'y n + 
z n {(pn + c n Zn) + z n d in terms of adjustment error 6 and adjusted parameter vector 9, and developed 
the following steps: 

i) Introduce the constant \/b mn in the parameterization, see (22). 

ii) Introduce the inequality \]b mn < y\b\, see (27). 

Hi) Apply the Young's inequality to obtain b, see (29). 

Recall that the value ofb mn is not required to be known. 

Substituting (30) into (20), we obtain: 

z„z„ < -c„zl + {3/A)C bvz + bz n (u + jJ— sgn(b)z„(<p T 0) 2 J 

-VK^\z„\y T 6 
we choose the following control law: 

u = --^s S n(b)z n (cp T §) 2 (32) 

3^-bvz 

where q>, z n are defined in (23), (16), respectively. Substituting (32) into (31), we obtain: 

z„z„ < -c n z\ + (3/4)C te2 - \ / % m \z n \f T 6 (33) 

To handle the effect of the constant (3/4)C/,; J2 , we formulate the following Lyapunov-like 
function: 

((l/2)(VV z --i/€ b ^) 2 UV z >C bvz 
V z = \ (34) 

| otherwise 

V 2 = (l/2)( 2 2 + ...+ 2 2) (35) 

where Cj, vz is defined in (28). We need the following properties: 
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Proposition 5.1. The function V z defined in (34) has the following properties: 

i)V z > (36) 

ii)V z <3C hvz +3V z (37) 

iii)V z , (dV z I dV z ) are continuous with respect to V z (38) 

Proof. From (34) it follows that V z > OVt > t , the property i of proposition 5.1. In addition, 
from (34) it follows that 

V 2 <(V2^+v/Q^) 2 (39) 

Applying the Young's inequality (cf. (Royden, 1988) pp. 123), we obtain: 

V z = C bvz + 2^/C^^W, + 2V Z < 3C bvz + 3V Z (40) 

This completes the proof of property ii. From (42) it follows that dV z /dV z = if V z = C bvz 
and that dV z /dV z is continuous. From (34) it follows that V z is continuous. This completes the 
proof of property iii of proposition 5.1. 



□ 



Differentiating (34) with respect to time, we obtain: 

dV z dV z 



dt w, v * (41) 



d% = r (1/2)(WV 2 )(VV Z - Vc bvz ) if v z > c bvz 

dV 7 } otherwise 



(42) 



To compute V z , we differentiate (35) with respect to time: V z = ziZi + ■ ■ ■ + z n z n . Introducing 
(11) and (15), we obtain 

V z = -c x z\ + ZJZ2 - c 2 z\ + z 2 zt, H h z„z„ (43) 

substituting (33), we obtain: 

V z < -ciz\ + zxz 2 - c 2 z\ + z 2 z 3 H c n z\ 

+ (3/4)C te2 - yt|z„|y T 6 
Provided that C\ > 2, c 2 > 2, ■ ■ ■ , c n > 2 and completing the squares yields: 
-C\z\ + zjz 2 - c 2 z\ + z 2 z 3 ■ ■ ■ - c„Z* 



(44) 



< 



z{ - (3/4)^ + (3/4)zf, < - (3/2) V z 



substituting into (44), we obtain: 

V z < -(3/2)V z + (3/4)C bDZ - VK;,\z„\cp T e (45) 
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Since dV z /dV z is non-negative, we can multiply it into (45) without changing the direction of 
the inequality: 

wj> * -t 3/2 ^w z + w^w z - ^™m<p t °w z (46) 

Substituting into (41), we obtain: 

dVy dV, dV 7 , T ~dV 7 

-£ < -0/2)V,^ + (3/4)C to2 ^ - ^KnW^O^ (47) 

we choose the update law so as to reject the effect of the term involving the adjustment error 
9: 

A dV-r 

e = Iq>\z n \^ (48) 

where T is a diagonal matrix whose diagonal elements are positive constants defined by the 
user, whereas q>, z n , dV z /dV z are defined in (23), (16), (42), respectively. 

Remark. So far, we have developed the controller, which involves the control law (32) and the update 
law (48). Other parameters necessary for its implementation are: V z defined in (35); Z\, Z%, ■ ■ ■ , Zn 
defined in (9), (12), ■ ■ ■ , (16), respectively; C\, vz defined in (28). Recall that c\ > 2, ■ ■ ■ , c n > 2 are 
user-defined positive constants. 

Remark. The control and update laws stated in (32) and (48) have the following features: 

i) The control law uses the adjusted parameter vector 8, so that it does not rely on upper or lower bounds 
of the plant coefficients, i.e. \i\, ■ ■ ■ , p.;, u&, b mn , and excessive control effort is also avoided. 

ii) Additional states are not required to handle the unknown varying behavior of the control gain, what 
is an important benefit with respect to closely related schemes that use the Nussbaum gain method. 

Hi) The control and update laws do not involve discontinuous signals. In fact, the vectorial field of the 
closed loop system is Lipschitz continuous, so that trajectory unicity is preserved. 

6. Boundedness analysis 

In this section we prove that the closed loop signals Z\, ■ ■ ■ , z n , 6, u are bounded if the 
developed controller is applied. 

Theorem 6.1. Boundedness of the closed loop signals. Consider the plant (1) subject to 
assumptions Ai to Av; the signals Z\, • • • , z„ defined in (9), (12) and (16); the signals <p, V z , dV z /dV z , 
Cb vz defined in (23), (35), (42) and (28), respectively. If the controller (32), (48) is applied, then the 
signals Z\, ■ ■ ■ , z n , 8, and u remain bounded and \e\ is upper bounded as follows: 

\e\ < V2 (VC^+ \j2V{x{t ))^j (49) 

Proof. We choose the following Lyapunov-like function: 

V(x(t)) = V z + V 9 (50) 
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Vg = {l/2)^K„e T I- l 8 (51) 

X(t) = [z v .... Z„, 9 T ] T (52) 

where V z is defined in (34) and 8 in (26). The time derivative of Vg is: 

v e = (i/2) v /fo^(e T r- 1 + e T r- 1 6i) (53) 

Since T is diagonal, then T _1 is diagonal, (T _1 ) T = T _1 , T r _1 = T r _1 0. In view of this 
and the update law (48), we have: 

Vg = y/bmnFT-H = Vb^,e T flZn^ (54) 

Differentiating (50) with respect to time, we obtain: V = V z + Vg. Substituting equations (47) 
and (54), we obtain: 



V<-(3/2)V z ^ + (3/4)C tez ^ 

_ _33]4 (Vz , ]4 _ c bm \ 
~ 2 3V 2 ^ 2 + 2 2 ) 



(55) 



From (42) if follows that 



dV z /dV z = for V z < C bvz (56) 

dV z /dV z > for V z > C bvz . (57) 



In view of this and (55), we obtain: 



. 3dV z (V z 



3 dV z (Vz 
2 dV z \ 2 



^<o = -o5^hf )^v z <c bvz 



^ V ^-Tw z v ^° (58) 

Thus, V + 0V < 0. Using the Lemma in (Slotine & Li, 1991) pp. 91, we obtain: 

V(*(*)) < V(x(t o ))exp(-0t) = V(x(t )) (59) 



where 



V(x(t )) = V Z0 + Vg (60) 

f (l/2)(vT^- VQ^) 2 if v 20 > C tez 

V z „ = { (61) 
[ otherwise 

V zo = (l/2)(2 1 (f ) 2 + ---+2„(f ) 2 ) (62) 

V fl0 = (1/2) VbmniHto) - e^T' 1 ^) - 6) (63) 
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Since V{x(t)) > 0, we have: < V(x(t)) < V(x(t )). Introducing the definition (50), we 
obtain 

V z + V <V(x(t o )) (64) 

=*• % < V(x(t )), V e < V(x(t )) (65) 

Thus, it follows from (51) that 6 £ Loo, and consequently 6 £ Leo. The inequality V z < V(x(t )) 
implies that the tracking error e is bounded, as we show hereafter. We begin by solving (34) 
forV z : 

V z = (^C^+VzVzf if V z > (66) 

V z < C bvz HV z = (67) 

Using the inequality V z < V (x(t )), we obtain: 



Vz < ( VC bvz + ^2V(x(t )) ) if V z > (68) 

\ -_ < C bvz < ( ^C bvz + j2V(x(t )) ) if V z = (69) 



combining both inequalities, we obtain: 

Vz < (Vc^ + y/2V(X(t ))\ (70) 



Introducing the definition (35), we obtain: 



lz\ + ■ • • + Z 2 < V2 f VQ^+ V2V(*(f ))J (71) 

so that Z\ 6 Loo, ■ ■ ■ , z„ G Loo- Since e 2 = z\ < z\ + ■ ■ ■ + z\, we obtain: 

\t\ < s/2 (VCbVz + \/2V(x(t ))) (72) 



which indicates the upper bound for the tracking error e. 

Remark. Notice that this upper bound does not involve integral terms, what is an important advantage 
with respect to the Nussbaum Gain method, see (Su et ah, 2009), (Feng, Hong, Chen & Su, 2008), (Feng 
et al., 2006), (Ge & Wang, 2003), (Feng et al., 2007), (Feng, Su & Hong, 2008), (Ren et al., 2008). 

Now, we proceed to show the boundedness of u. From (9), (12), (16) it follows that X\ 6 Loo, 
%2 £ i-oo, ■ ■ ■ , x n g Loo- Therefore, 7,, € Loo- It follows from (23) that y 6 L M . From (32) it 
follows that u 6 Loo- This completes the proof. 

□ 
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7. Convergence analysis 

In this section we prove that if the developed controller is applied, then the signal Z\ converges 
asymptotically to 2 , where Q 2 = {z\ '.\z\\ < C be }- 

Theorem 7.1. Convergence of the tracking error. Consider the plant (1) subject to assumptions Ai 
to Av; the signals z\, ■ ■ ■ , z„ defined in (9), (12) and (16); the signals f, V z , dV z /dV Z/ Q, ra defined in 
(23), (35), (42) and (28), respectively. If the controller (32), (48) is applied, then the signal Z\ converges 
asymptotically to fi z , where fi z = {z\ : |zi| < C be }. 

Proof. In view of (42), equation (58) can be rewritten as: 
V < 
fd = 



fd<0 




(73) 


' (3/8)WV~ z )(y/V z - 



- vQroz) if V z > C bvz 

otherwise 


(74) 



The derivative dfy/dVz is not continuous, as it involves an abrupt change at V = C bvz . Thus, 
the Barbalat's Lemma can not be applied on fy. To remedy that, we shall express (73) in terms 
of a function with continuous derivative: 

V < -ft <-fg<0 (75) 

_ f (3/8) (v 7 ^- v/C^) 2 if V z > C bvz 
f S-{ otherwise { ' 

Arranging and integrating (75), we obtain: 

fg<~V 



1 f s dT< V(x(t ))-V(x(t)) 



L 

V(x(t))+ f'f g dT<V(x(t )) (77) 

J to 

Thus, fg £ L\. We have to prove that fg e L^, f„ £ L^ to apply the Barbalat's Lemma. Since 
Vz 6 Loo, it follows from (76) that fg G Leo. Differentiating (76) with respect to time, we obtain: 

t* = wk v * _ _ _ (78) 

% = f (3/8)(l/v / V 2 )( v / V z - VQ^) if V z > C bvz 

dV z \ otherwise V ' 

Notice that df g /dV z is continuous. Since V z 6 Loo, then dfg/dV z 6 Loo- Since Z-y G Loo, • • • , 
z n G Loo, it follows from (11), (15) that i\ 6 Loo, ■ ■ ■ , i n —l € Loo- Since u 6 Loo, it follows from 
(17) that z n e Loo- Therefore, from (43) it follows that V z G Loo- 

So far we have proved that dfg/dV z € Loo and V z 6 Loo, so that it follows from (78) that 
fg £ Loo- In view of fg £ Loo, fg £ Loo, application of Barbalat's Lemma (cf. (Ioannou & 
Sun, 1996) pp. 76), then indicates that fg converges asymptotically to zero. Hence, from (76) 
it follows that V z converges to Ci vz , where Cl vz = {V z : V z < C bvz }. From the definition (35), 
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it follows that Z\ converges asymptotically to Q z , where fi z = {zi : \l\\ < ^/2C], VZ }. Since 
Cbvz = (1/2)C? / it follows that fi 2 = {zi : \z\\ < Q, e }. This completes the proof. □ 



8. Simulation example 

Consider the following case of the plant (1): 

y = y 2 a + bu + d 

72 = [y, y} T , a = [a x , a 2 ] T 

fll = -2(l + 0.1sin(2(7r/8)f)), a 2 = 
b = 2 (1 + 0.1sin((27r/ll)f)) + 0.6|y| 
d = -0.2 (1 + 0.1sin((27r/7)f)) y 



-l(l + 0.1sin(2(7r/5)f)) 



(80) 
(81) 

(82) 



The aim is that y converges towards y^, with a threshold of 0.1. In figure 1 we present a 
simulation block diagram for the example. 





ys.ftf.&i 








) 










Update 
Law 










A 
























Reference 
Model 


1 — - 


Control 

Law 


u 


Nonlinear 
System 


y,y 




(—> 

























Fig. 1. Simulation block diagram. 

The properties Ai, Aiv, Av of section 2 are analyzed at the following. From (82) it follows that 
a.\,d,b are bounded as: 



|fli| < 2(1.1) = 2.2, |a 2 | < 1(1.1) = 1.1 
|d|< 0.2(1.1) |y|=0.22|y| 
\b\ > 2(0.9) = 1.8 > 

Hence, the upper bounds of fli , a 2 , d, and the lower bound of b, are: 

[«i| < H = 2.2, \a 2 \ < }i 2 = 1.1, \b\ > b mn = 1, 
M < Fd/d/ f'rf = 0-22, f d = \y\ 



(83) 
(84) 
(85) 



(86) 
(87) 



where /^ is not constant and known, whereas b mn , \i\, fi 2 , fi^, b mn are positive, constant and 
unknown to the controller. From (86), (87) it follows that assumptions Ai, Aiv, Av of section 2 
are satisfied. 

The procedure of section 4 is followed in order to establish the terms involved in the control 
and update laws, mentioned in remark 5. Eq. (80) can be rewritten as: 



X\ = x 2 

x 2 = y 2 a + bu + d 

x\= y, x 2 = y, n = 2 



(88) 
(89) 
(90) 
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since n = 2, the state transformation based on the backstepping procedure involves the steps 
0, 1, 2. 

Step 0. Let 

z l =e = y-y d = x 1 -y d (91) 

as in (9). 

Step 1. Differentiating (91) with respect to time and arranging, yields: 

zj = -c lZl + z 2 (92) 

z 2 = x 2 + c x z\ - y d (93) 

as in (11), (12). 

Step 2. Since n = 2, the second step is the last one. Differentiating (93) with respect to time, 
using (89) and arranging, yields: 

z 2 = x 2 + c x Z\ - y d 

= yja + bu + d + c\i\ - y d 

= j 2 T a + bu + d + c 1 (x 2 + cp 1 )-y d (94) 

using the definitions (91), (93), yields: 

z 2 = yja + bu + d+ cp 2 (95) 

9>2 = ci(z 2 -ciZi)-y d (96) 

notice that the form of (95), (96) is that of (17), (18), respectively. This completes the state 
transformation based on the backstepping procedure. 

The parameters defined above can be summarized as: 

zi=y-y d (97) 

z 2 = x 2 + c 1 z 1 - y d (98) 

*i = y, x 2 = Vd (99) 

n = -n (ioo) 

n = ci{z 2 -c 1 zi)-y d (101) 

According to remark 5, it remains to define <p, V z . From (81), definition (23) and n = 2, it 
follows that 

?= [|72[1]I/ lT2[2]|/ fd, !<P2 + c 2 z 2 |] = [\y\, \y\, fd, \f2 + c 2 z 2 \] T (102) 

From (35) and n = 2 it follows that 

V z = (l/2)(zl + z 2 2 ) (103) 
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Expressions (97) to (103) allow to define the control and update law. From (32), (48), (82) and 
n = 2 it follows that 



sgn(fo) 
u = — 



= +1 
1 

dVz 



z 2 (cp , e) 2 



Tf\z 2 



(104) 
(105) 

(106) 



the main parameters needed to compute u and 6 are: <p (102), q> 2 (101), Cbvz (28), z 2 (98), Z\ 
(97), dV z /dV z (42), V 2 (103). In addition, T is a diagonal matrix whose diagonal elements are 
positive constants defined by the user. 




Fig. 2. Example 1, upper: output y (continuous line), desired output y^ (dash-dot line); 
middle: tracking error e; lower: control input u. 

Since the aim is that y converges towards y^, with a threshold of 0.1, we set Q, c = 
We use the reference model (5) with yd{t ) = y{t ), Vdi^o) = 0, a, % i = 1, a mp 
We use the following parameter values for the control and update laws: C\ = 2, c 2 
r = diagjl, 1, 1, 1}. 



0.1. 
= 1. 

= 2, 



The results are shown in figures 2 and 3. We have choosen y^(fo) ~ 3/(^0 ) m order to obtain a 
rapid convergence of y towards y^. Figure 2 shows that. /') the tracking error e converges 
asymptotically towards fi e = {e : |e| < 0.1}. ii) The output y converges towards y^ 
with threshold 0.1 without large transient differences. Figure 3 shows that §i,...,9^ are not 
decreasing with respect to time. This occurs because 8 is non-negative. The procedure for the 
sample plant (80) is simpler in comparison with adapive controllers that use the Nussbaum 
gain method. 
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1. Introduction 

Piezoelectric actuators are the most suited actuation devices for high precision motion 
operations in the positioning tasks include miro/nano-positioning [1]. These actuators have 
unlimited motion resolution and posses some advantages such as ignorable friction, 
noiseless, zero backlash and easy maintenance, in comparison with the conventional 
actuated systems which are based on the sliding or revolute lower pairs [2] . 

Producing large forces, fast response and high efficiency are major advantages of 
piezoelectric actuators. But, it has some drawbacks such as hysteresis behavior, drift in time, 
temperature dependence and vibration effects. Molecular friction at sites of materials 
imperfections due to domain walls motion is the general cause of hysteresis in piezoelectric 
materials [3]. The hysteresis is a major nonlinearity for piezo-actuators and often limits 
system performance via undesirable oscillations or instability. Therefore, it is difficult to 
obtain an accurate trajectory tracking control. Numerous mathematical methods have been 
proposed to analyze the hysteresis behavior of piezoelectric actuators. These studies may be 
categorized in asymmetrical and symmetrical methods. The asymmetrical types of hysteretic 
models include polynomial model [4], Preisach's model [5], neural network model [6] and 
Karasnoselskii and Pokrovskii [7]. The symmetrical types of hysteretic model include 
Duhem model [8], Bouc-Wen model [9] and Lugre model [9]. 

The asymmetrical methods establish the nonlinear relations between the input and output 
based on the measured input/ output data sets. Because superposition of a basic hysteresis 
operator is a fundamental principle in these models, they are also called to operator based 
model. Although, an operator based model may give a good match with experimental data, 
the dynamics of the piezoelectric material is not formulated in these modeling methods and 
model parameter identification and implementation is more difficult in this case. The 
symmetrical methods employ nonlinear differential equations in order to describe 
hysteresis. In this case, the dynamics of the piezoelectric materials are described but the non- 
symmetric hysteresis is not modeled. However these models are more tractable for control 
design. In order to include the hysteresis effect and compensating its effect, Lugre model 
[10] is analyzed and studied in this chapter. 
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There are several control methods to overcome the above mentioned errors and increase the 
tracking control precision of the piezoelectric actuators. Some of these methods are PI and 
PID controller, fuzzy controller [11], adaptive RFNN [12], feed-forward model reference 
control method [13], adaptive hysteresis inverse cascade with the plant [14], reinforcement 
discrete neuro-adaptive controller [15], adaptive wavelet neural network controller [16], 
nonlinear observer-based sliding-mode controller [17] an adaptive backstepping controller 
[10, 18], robust motion tracking controller based on sliding-mode theory [19] and continuous 
time controller based on SMC and disturbance observer [20]. In some of these works, a 
complex inverse hysteresis model has been adopted to overcome the nonlinear hysteresis 
effect. Also, in the methods based on the neural network approach, to ensure the error is 
bounded, it is assumed that the system states must be inside a compact set. Moreover, 
robustness against parameters uncertainties and external disturbances is the other problem 
encounter the control methods presented for piezoelectric actuator. 

In this chapter, a robust motion tracking control strategy in combination with the hysteresis 
force observer is designed and investigated for the piezoelectric actuator. The presented 
controller is robust against the unknown or uncertain system parameters and can estimate 
the hysteresis force with its estimation property. This control strategy is established based 
on the lumped parameter dynamic model. Using Lyapanouv stability analysis, the stability 
analysis of the overall control and observer system is performed. Furthermore, the validity 
and effectiveness of the designed methodology is investigated by numerical analysis and its 
results obtained are compared with those of [19] . 

Section 2 contains the model explanation of piezo-positioning mechanism. A hysteresis 
model for piezoelectric systems based on Lugre model is discussed in section 3. Design 
procedure of the developed controller is presented in section 4. Simulation analysis and 
results obtained are presented in section 5. Finally, section 6 includes the conclusion of this 
chapter. 

2. Model of piezo-positioning mechanism with hysteresis 

The lumped parameter dynamic model of the piezo-driven mechanism is written as [10, 16]: 

M'x + Dx + F H +F L =u (1) 

Where M is the mass of the controlled piezo-positioning mechanism, D is the linear friction 
coefficient of the piezo-driven system, F L denotes the external load, F H is the hysteresis 
friction force function, x(t),x(t), and x(t) denote the piezoelectric displacement, velocity and 
acceleration, respectively and u is the applied voltage to the piezo-positioning mechanism. 

A block diagram of the model (1) is depicted in Fig.l. Note that the simulation program 
used in this chapter is C++. 

Noted that: 

M0| * ^ H0 , \P H (t)\ < SF m , |F H (f)| < SF H2 „ \F L (t)\ < SF L0 \F L (t)\ < 8F L1 , \P L (t)\ < SF L2 (2) 
Where SF Hi e 9? ,i = 0,1,2 and SF U e 5R ,i = 0,1,2 denote the known upper bounds. 
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Fig. 1. Model of piezo-positioning mechanism with hysteresis [10]. 

3. Lugre hysteresis model for piezo-positioning system 

This chapter uses Lugre model to describe nonlinear hysteretic curve of piezoelectric 
systems. The mathematical equation is as follows [10,16]: 



F H = o n z -a - -, zljcl+foi + a 7 )x 

l g{x) M l 2 > 



(3) 



z=x — — — z 
8(x) 



(4) 



Where z interpreted as the contact force applied voltage average bristle coefficient, 
<T ,<T 1 ,a" 2 are positive constants and generally can be equivalently interpreted a bristle 
stiffness, damping and viscous-damping coefficients, respectively. Moreover, the function 
g(x) denotes the Stribeck effect curve described by: 



°og(x) = fc+(fs-fcV msf 



(5) 



Where f c is the coloumb friction level, f s is the level of stiction force and x s is the Stribeck 
velocity. The function g(x) is positive and constant [16]. As depicted in [10], the following 
lemma is held. 

Lemma 1. Consider nonlinear dynamics system of (4). For any piecewise continuous signal 
x and x , the output z(f) is bounded. 

Substituting the (3) in (1) and arranging the expression, the following dynamics equation is 
obtained for the piezo-positioning mechanism with Lugre hysteresis model: 
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4. Design controller 



In this section we design a robust control strategy for the system of (1), to asymptotically 
estimate the hysteresis parameter of equation (3). To facilitate the design process, we 
assumed that the hysteresis force is dependent only on the time and F H and its first two 
time derivatives remain bounded for all times. Using (6) we can rewrite equation (1) as: 

.. u F D . 

x = x (7) 

M M M 

z x 
Where F = (a z + F L )-a 1 — — + (<r 1 +a 2 )x is a combination of an unknown friction 

hysteresis force function and external load, which must be estimated. 

The dynamics of a piezo-positioning system can be represented by the following equation: 

Mx + Dx + F-m + (AMx + ADx) = (8) 

Where M and D are the nominal parameter values of the mass and linear friction coefficient 
and AM and AD are the parametric errors between real value and nominal value of the 
uncertain parameters of the system and modeled as: 

|AM| < SM (9) 

|AD| < SD (10) 

Where SM and SD are the bounds of system parameters. The position tracking error 
signal is defined as: 

e = x d -x (11) 

Where x d s SR denotes the desired position trajectory. The desired piezo-position trajectory 
and its first three time derivatives are assumed to be constrained by the following: 

N < £d0'|*d| < £d\>\*d\ < Cd2>\ X 'd\ < Cd3 (12) 

Where the Cdi s denote known positive constants. Also, the filtered tracking error signal 
r(t) e 3? 1 is defined as: 

r = e + ae (13) 

Where a e SR is a positive constant control gain. 

The system dynamics of (8) are rewritten in terms of the filtered tracking error signal r(t) as 
follows: 
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v d ' M M M M 

We define % as: 

% = ^L(AMx + ADx)<^L(SM\x\ + SD\x\) = ^Lfi 1 (x) (15) 

Another filtered tracking error signal s(i) e 9? is defined as: 

s = r + fir (16) 

Where /? is a constant positive parameter. 

Based on the developed error system and ensuring the stability analysis, the following 
control input signal for the system (14) is designed: 

u = M(x d + ae) + Dx + F + /7 1 (x)sgn(s) (17) 

Where F e SR represents the estimation of F and is obtained by the following observer: 

F = -(k 1 +j3)F + k 1 j3r + psgn(r) (18) 

Where sgn(.) is the standard signum function, and k lr p e s Jl are positive constants. 

Noted that the equation (18) for F is a stable linear system with the disturbance 
term k-^fir + /?sgn(r) . 

To facilitate the dynamic system stability analysis, the auxiliary disturbance signal 
t](t) e 91 is defined by: 

r 1 = F + (k 1 +p)F (19) 

Due to the boundness of F and its first two time derivatives, it is understand that 
7(t),i7(t)eL.,. 

Based on the stability analysis, the constant p should be chosen to satisfy the following 
inequality: 

p>\rj(t)\ + ±\v(t)\ (20) 

By substituting (17) in (14), and simplify the obtained expression, the dynamics of f is 
achieved as: 

r=F-F + i(*i(*)-A(*)sgn<s)) (21) 

Where Xi{ x ) = AM* + ADx . 
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Now, using (21), the time derivative of (16) is obtained as: 

s = F-F + i(ZiW-AWsgn(s)) + /7(F-F + i(z 1 (x)-AWsgn(s))) (22) 

Substituting (18) and (19) in (22), gives: 

s = jj - k lS - psgn(r) + -j^(Xi(x) - A(*)sgn(s)) + — ^"UiM - A(*)sgn(s)) (23) 

Remark 1: If s(f) s L^ , then r(t),r(t) eL x and if s(f) is asymptotically regulated, then 
r(t),f(t) are also, asymptotically regulated. 

5. Stability analysis 

Theorem 1. For the dynamics of (7), the designed controller of (17) and (18) guarantees the 
global asymptotetic piezo-driven position tracking in the sense that: 

lime(i) = (24) 

f->cc 

And global asymptotetic estimation of hysteresis in the sense that: 

lim[F-F] = (25) 

With the constant p satisfies the condition of (20). 

Proof 

A non-negative, scalar function V(t) s s Jt is defined as: 

V=|s 2 (26) 

Tacking time derivative of (26) gives: 

y = -lc 1 s 2 + (r + /7r)( 7 -psgn(r)) + -L(i' 1 (x)-AWsgn(s)) + ^i r is( Zl (x)-A(x)sgn(s)) (27) 

M M 

From (15) it is clear that: 

i(Z 1 (x)-AWsgn(s)) + ^(^ 1 (x)-AWsgn(s))<0 (28) 

Where t (x) = SM\x\ + SD\x\ . 
Hence: 

V <-k r s 2 +(f + /3r)(r/-psgn(r)) (29) 

After integrating both sides of (29), the following inequality is obtained: 
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j*??(g) 



da 



- p Ida 



+ Co (30) 



V(t) -V(t )< -kX s 2 {a)da + [\r(t)\ \tft)\ - p\r(t)\] + f ' p\r{*)\ L\ + \ 

Where £" G s -^ ^ s a positive constant, defined by: 

fo = k(*o)||'/(*o)| + /'K*o)| ( 31 ) 

After applying the (20) to the bracketed term of (30), V(t) can be upper bounded as follows: 

V(t)<V(t )-kJ t S 2 (a)da + C (32) 

It is deduced from (32) and (26) that V(t ) e L M and s(t ) e L^ , respectively. Finally utilizing 
remark 1, r(i),f (f),e(i),e(i),s(f) e L„ . 



Equation (32) rearrange as: 



U f sV)da<V(f )-V(f) + ^o 

Jf 



(33) 



From the fact that V(t) is non-negative and equation (33), it can be deduced that s(t) s L 2 . 
Because of s(f)eL oo ni- 2 and s(f) s L K , we can use the Barbalat's Lemma [21] to 
summarized that: 



Therefore: 



lims(f) = 



limr(f),r(r),e(r),e(r) = 

£-»co 



(34) 



(35) 



6. Simulation results 

In this section, simulation results are presented to investigate the performance of the 
presented method for piezoelectric actuator. First, we test the controller response with the 
nominal value of piezoelectric actuator, when F L =0.3sin(f) and initial values are: 
x(0) = 0,x(0) = 0,z(0) = . The objective of the positioning is to drive the displacement signal 
x to track the reference trajectory which is shown in Fig.2. 

The parameters of piezo-positioning mechanism are given in Table 1. 



<r = 50000N / m 




a j = 0.4Ns / m 


o- 1 =V5xl0 4 Ns/m 


f c =lN 


f s =1.5N 


x s = 0.001m / s 


M = lkg 


D = 0.0015Ns / m 


F L =0.3sin(f) 



Table 1. Piezoelectric Parameters 
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Fig. 2. Desired and Actual Piezoelectric displacement. 
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As depicted from the results, especially Fig.4, the positioning mechanism has a good 
accuracy and the error between displacement and reference signal at least 4 orders is smaller 
than the actual signal which means that the error is about %0.01. 
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Fig. 4. Error between actual & desired displacement. 

For the purpose of study the behavior of the controller in the presence of parameter 
uncertainties, by considering theAM=0.1M, AD = 0.1D and AF L =0.4and using 
x(0) = 15(//m),x(0) = 0,z(0) = 0as the initial values, we have performed another test. The 
result of this test is compared with the result of the method presented in [19]. In the 
simulation we consider the control gain a is 2000. 

As it can be seen from Fig. 5 and Fig. 7, the presented method has an acceptable response 
and the positioning error is in about %0.25. Also, the hysteresis and disturbance voltage 
and its estimation are shown in Fig. 9. As it can be seen from this result, the hysteresis 
identifier can estimate the hysteresis voltage precisely. For the purpose of comparison the 
accuracy of the proposed method with the recently proposed method in [19], we simulate 
the piezo-positioning mechanism with the method of [19]. The displacement error of [19] 
is shown in Fig. 10. Noted that in the simulation of [19] we use these control gains: 
k =W 5 ,k v =9000,fc s =50 and a =10. Comparison of Fig. 7 and Fig.10 depicts that the 
error in the presented method is smaller than the method of [19]. Also, it is clear that in 
addition of smaller number of control gains, the gain of the controller in the presented 
method in comparison with [19] is very smaller. Experimental implementation of the high 
gain needs more complexity and also it may be generate noise. Although method of [19] is 
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robust against disturbances and uncertainties and in comparison with other methods has 
higher precision, it needs high value gains to perform this task. While the proposed 
method can perform these tasks with lower cost and also it has the capability of hysteresis 
estimation. 

For the comparison between the presented method and the method of [19] another 
simulation test is carried out. In this case AM = 0.2M , AD = 0.2D and AF L = 0.4 . The result is 
shown in Fig. 11. 

As it can be depicted the presented method has the better response and using the proposed 
controller the error is decreased. Note that in this case the reference is similar to the 
Fig.5. 

The last simulation is devoted to the reference signal of Fig.2 by considering AM = 0.2M , 
AD = 0.2D and AF L = 0.4 . Also in this case the error of the proposed controller is less than 
the method of [19]. 
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Fig. 11. Comparison between the displacement error of the presented controller and the 
controller proposed in [19]. 
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Fig. 12. Comparison between the displacement error of the presented controller and the 
controller proposed in [19]. 

7. Conclusion 

A robust tracking motion controller is designed to control a positioning of piezoelectric 
actuator system with hysteresis phenomenon. The Lugre hysteresis model is used to model 
the nonlinearities in the system under study. Using the presented controller, we can 
estimate the nonlinear hysteresis and disturbances imposed to the piezoelectric actuator and 
compensate their effects. The performance and efficiency of the designed controller for a 
positioning system is compared with recently proposed robust method. The results obtained 
depict the validity and performance of the presented approach. 
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1. Introduction 

Control allocation is the process of mapping virtual control inputs (such as torque and force) 
into actual actuator deflections in the design of control systems (Benosman et al v 2009; 
Bodson, 2002; Buffington et al., 1998; Liao et al., 2007; 2010). Essentially, it is considered as 
a constrained optimization problem as one usually wants to fully utilize all actuators in order 
to minimize power consumption, drag and other costs related to the use of control, subject to 
constraints such as actuator position and rate limits. In the design of control allocation, full 
state information is required. However, in practice, states may not be measurable. Hence, 
estimation of these unmeasurable states becomes inevitable. 

The unmeasurable states are generally estimated based on available measurements and 
the knowledge of the physical system. For linear systems, the property of observability 
guarantees the existence of an observer. Luenberger or Kalman observers are known to give 
a systematic solution (Luenberger, 1964). In the case of nonlinear systems, observability in 
general depends on the input of the system. In other words, observability of a nonlinear 
system does not exclude the existence of inputs for which two distinct initial states generate 
identical measured outputs. Hence, in general, observer gains can be expected to depend on 
the applied input (Nijmeijer & Fossen, 1999). This makes the design of a nonlinear observer 
for a general nonlinear system a challenging problem. Although various results have been 
proposed over the past decades (Ahmed-Ali & Lamnabhi-Lagarrigue, 1999; Alamir, 1999; 
Besancon, 2007; Besancon & Ticlea, 2007; Bestle & Zeitz, 1983; Bornard & Hammouri, 1991; 
Gauthier & Kupka, 1994; Krener & Isidori, 1983; Krener & Respondek, 1985; Michalska & 
Mayne, 1995; Nijmeijer & Fossen, 1999; Teel & Praly, 1994; Tsinias, 1989; 1990; Zimmer, 1994), 
none of them can claim to provide a general solution with the same convergence properties as 
in the linear case. 

Over the past decades, a variety of methods have been developed for constructing nonlinear 
observers for nonlinear systems (Ahmed-Ali & Lamnabhi-Lagarrigue, 1999; Alamir, 1999; 
Besancon, 2007; Besancon & Ticlea, 2007; Bestle & Zeitz, 1983; Bornard & Hammouri, 1991; 
Gauthier & Kupka, 1994; Krener & Isidori, 1983; Krener & Respondek, 1985; Michalska & 
Mayne, 1995; Nijmeijer & Fossen, 1999; Teel & Praly, 1994; Tsinias, 1989; 1990; Zimmer, 
1994). They may be classified into optimization-based methods (Alamir, 1999; Michalska 
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& Mayne, 1995; Zimmer, 1994) and feedback-based methods (Bestle & Zeitz, 1983; Bornard 
& Hammouri, 1991; Gauthier & Kupka, 1994; Krener & Isidori, 1983; Krener & Respondek, 
1985; Teel & Praly, 1994; Tsinias, 1989; 1990). Optimization-based methods obtain an estimate 
x(t) of the state x(t) by searching for the best estimate x(0) of x(0) (which can explain 
the evolution i/(t) over [0, (]) and integrating the deterministic nonlinear system from x (0) 
and under m(t). These methods take advantage of their systematic formulation, but suffer 
from usual drawbacks of nonlinear optimization (like computation burden, local minima, 
and so on). Feedback-based methods can correct on-line the estimation x(t) from the 
error between the measurement output and the estimated output. These methods include 
linearization methods (Bestle & Zeitz, 1983; Krener & Isidori, 1983; Krener & Respondek, 
1985), Lyapunov-based approaches (Tsinias, 1989; 1990), sliding mode observer approaches 
(Ahmed-Ali & Lamnabhi-Lagarrigue, 1999) and high gain observer approaches (Bornard & 
Hammouri, 1991; Gauthier & Kupka, 1994; Teel & Praly, 1994), and so on. Among them, 
linearization methods (Krener & Isidori, 1983) transform nonlinear systems into linear systems 
by change of state variables and output injection. It is applicable to a special class of nonlinear 
systems. Sliding mode observer approaches (Ahmed-Ali & Lamnabhi-Lagarrigue, 1999) is 
to force the estimation error to join a stabilizing variety. The difficulty is to find a variety 
attainable and having this property. High gain observer approaches (Besancon, 2007) use the 
uniform observability and weight a gain based on the linear part so as to make the linear 
dynamics of the observer error to dominate the nonlinear one. Due to the requirement of the 
uniform observability, these approaches can only be applied to a class of nonlinear systems 
with special structure. Interestingly, Lyapunov-based approaches (Tsinias, 1989; 1990) provide 
a general sufficient Lyapunov condition for the observer design of a general class of nonlinear 
systems and the proposed observer is a direct extension of Luenberger observer in linear case. 

In this chapter, we extend the control allocation approach developed in (Benosman et al., 2009; 
Liao et al., 2007; 2010) from state feedback to output feedback and adopt the Lyapunov-type 
observer for a general class of nonlinear systems in (Tsinias, 1989; 1990) to estimate the 
unmeasured states. Sufficient Lyapunov-like conditions in the form of the dynamic update 
law are proposed for the control allocation design via output feedback. The proposed 
approach ensures that the estimation error and its rate converge exponentially to zero as 
t — y +oo and the closed-loop system exponentially converges to the stable reference model 
as t — v +oo. The advantage of the proposed approach is that it is applicable to a wide class 
of nonlinear systems with unmeasurable states, and it is computational efficiency as it is not 
necessary to optimize the control allocation problem exactly at each time instant. 

This chapter is organized as follows. In Section 2, the observer-based control allocation 
problem is formulated where the control allocation design is based on the estimated states 
which exponentially converge to the true states as t — >■ +oo. In Section 3, the main result of 
the observer-based control allocation design is presented in the form of dynamic update law. 
An illustrative example is given in Section 4, followed by some conclusions in Section 5. 

Throughout this chapter, given a real map /(v, w), (v, w) 6 R" x R m , D v /(vq, Wq) denotes 
its derivative with respect to v at the point (vq, Wq). For given real map h(y) with v G R", 
Dh(vo) denotes its derivative with respect to v at the point vq. In addition, || ■ || represent the 
induced 2-norm. 
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2. Problem formulation 

Consider the following nonlinear system: 

fx = /(x,u) 

\y = fc(x) W 

where x g X C R" is the state vector with X a open subset of R", y £ R is the measurement 
output vector, and u 6 R' H is the control input vector satisfying the constraints 

u £ fi = I u = [u\ ui ■ ■ ■ u,„] \u_j < Uj < fij, i = 1,2,' ■ ■ ,m\ (2) 

with u = [«i «2 ' ' ' Km] an d u = [«i «2 ■ ■ ■ fij»] being vectors of lower and upper 
control limits, respectively. 

We assume that the system (1) satisfies the following assumption: 

Assumption 1. The function /(x,u) is smooth and the output function h(x) is continuously 
differentiable. 

Since control allocation need full state information, the state estimation for the system (1) is 
required. 

Consider a dynamic observer of the following form 

x = /(x,u)-3>(x,u)[y-fc(x)] (3) 

Define the error e as 

e = x x (4) 

To estimate the state x, we wish to design the mapping ^(x, u) such that the trajectory of e 
with the dynamics 

e = /(x, u) - /(x, u) + *(«, u) [y -/.(*)] (5) 

exponentially converges to zero as t — ¥ +oo, uniformly on u 6 fi, for every x(0) subject to 
e(0) = x(0) — x(0) near zero. 

The aim is to design a nonlinear control allocation law based on the state observer (3) such that 
a reference model that represents a predefined dynamics of the closed-loop system is tracked 
subject to the control constraint u £ fi. 

Given that the predefined dynamics of the closed-loop system is described by the following 
asymptotically stable reference model 

x = A d x + B d r (6) 

where A^ G R" x ", B d £ ]R" x "r anc j me reference r G R" r satisfy the following assumption. 

Assumption 2. A d is Hurwitz, and r e E C R" r is continuously differentiable where E is an open 
subset defined by: for each r e E, there exist x €L X and u G fi such that the system (1) matches the 
reference system (6). 
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Since the state x is immeasurable, the control allocation design is then based on its estimate 
x. In other words, we have to first choose the mapping <5(x, u) in (3) such that the estimation 
error e exponentially converges to zero as t — > +oo, uniformly on u S CI, for every x(0) £ X 
subject to e(0) near zero; then minimize the cost function 

/(x,r,u) = -u T H 1 u+ -T r (x,r,u)H 2 T(x,r / u) (7) 

where Hi 6 ^rnxm anc j j_j 2 ^ R" x " are positive definite weighting matrices, and 

r(x,r,u) =/(x,u)- A d x - B d r (8) 

is the matching error between the actual dynamics and desired dynamics. Since power 

1 T 
consumption minimization introduced by the term -u Hi u is a secondary objective, we 

choose IIHjH -c || Hz || - 

Now the control allocation problem is formulated in terms of solving the following nonlinear 
static minimization problem: 



min / ( x, r, u ) subject to 
u£fl and x converges to x exponentially 



(9) 



Define 



with 



A(u) = [S( Ml ) S{u 2 ) ■ ■ ■ S(u m )} (10) 



S(Wj')=min((u,- m,) 3 , (m,-«,-) 3 ,0),;' = 1,2, ■ ■ ■ ,m (11) 

Then the constraint condition uGQis equivalent to 

A(u) = (12) 

Introduce the Lagrangian 

L(x,r,u,A)=/(x,r,u)+A(u)A (13) 

where A 6 R'" is a Lagrange multiplier. And assume that 



Assumption 3. There exists a constant 7j > such that — j > 7il« 



du 2 



The following lemma is immediate ((Wismer & Chattergy, 1978), p. 42). 

Lemma 1. If Assumptions 1 and 3 hold, the Lagrangian (13) achieves a local minimum if and only if 

— = and — = 0. 
dA du 
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Proof. Necessity: The necessary condition is obvious. Sufficiency: Since — = 0, we have 

dA 
A(u) = 0. In this case, the Lagrangian (13) is independent of the Lagrange multiplier A, which 

, , • • , 9L „ , d 2 L , d 2 L „ , , A 

achieves a local minimum if 3— = and t-^t > 0. As z—^ > is guaranteed by Assumption 

du du 2 du z 

3, thus, t-t- = and t-— = implies the local minimum. The proof is completed. □ 

dA du 

Remark 1. It should be noted that Assumption 3 is satisfied if all control inputs are within their 

limits (i.e., — - = A T (u) = 0) and the nonlinear system (1) is affine in control (i.e., /(x,u) = 
aA 

/l(x) + g(x)u). It is because, in this case, — ^ = Hj + g T (x)H2g(x) is positive definite matrix for 

du 
Hj > and H2 > 0. Furthermore, since the Lagrangian (13) is convex in this case, Lemma 1 holds 

for a global minimum. 

To solve the control allocation problem (9) with the state estimate x from the observer (3), we 
consider the following control Lyapunov-like function 

V(x, e, r, u, A) = V m (x, r, u, A) + i e T Pe (14) 

where P > is a known positive-definite matrix and 

1 



V m (x,r,u,A) 



dL\ T dL {dL\ T dL 



dxx) du + \d\J dA 



(15) 



Here the function V m is designed to attract (u, A) so as to minimize the Lagrangian (13). The 

1 T 
term - e Pe forms a standard Lyapunov-like function for observer estimation error e which 

is required to exponentially converge to zero as t — > +00. 

Following the observer design in (Tsinias, 1989), we define a neighborhood Q of zero with 
Q C X , a neighborhood W of X with {x- e : x e X, e e Q} C W, and a closed ball S of 
radius r > 0, centered at zero, such that S C Q. Then define the boundary of S as 3S. Figure 1 
illustrates the geometrical relationship of these defined sets. 

Let % denote the set of the continuously differentiable output mappings h(x) : X — > R such 
that for every mg £ Q and x g W, 

R(x,m )>0 (16) 

and 

kerR(x,m ) C kerD/?(x) (17) 

where 

R(x,m ) = [Dh(x)} T Dh(x + m ) + [Dh(x + m )] T Dh(x) (18) 

Remark 2. Obviously, every linear map y = Hx belongs to %. Furthermore, % contains a wide 
family of nonlinear mappings. 
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Fig. 1. Geometrical representation of sets 

We assume that 

Assumption 4. h(x) in the system (1) belongs to the set H, namely, h(x) G %. 

Further, we define 

A 



N 



je GR""|e r PD x /(x + mi ,u)e < -fc ||e|| 2 } 



(19) 



and assume that 



Assumption 5. There exist a positive definite matrix P g R"" x " x and a positive constant kg such 
that kerDh(x) C N holds for any (x,mi,u) G W x Q x n. 

Remark 3. Assumption 5 ensures that the estimation error system (5) is stable in the case ofh(x) = 
h(x) and x ^ x. In particular, for linear systems, the condition in Assumption 5 is equivalent to 
detectability. 



3. Main results 

Denote 



and define 



Let 



du 2 dAdu 3u 
d 2 L 



3u3A 







3A 



M 



A 



jveR"*| y = r||e|r 1 e, eGNns] 



7i(x,u) = max|r 2 (||P||||D x /(x + m 1/ u)|| +k ), m 1 e S, (x,u) G W x n| 
72 ( x ) = mm I ~v T R{%, mo)v, mp G S, v G dS — M, x G W 



(20) 

(21) 

(22) 
(23) 
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Theorem 1. Consider the system (1) with x £ X and u G CI. Suppose that Assumptions 
1-5 are satisfied. For a given asymptotically stable matrix A^ and a matrix Bj, given symmetric 

dL dL 



positive-definite matrices I" \ and I \,and a given positive constants to, ) 'or e(0) nearzero, \ —, — ,e 

\ a A au 

exponentially converges to zero as t — > +oo, and the dynamics of the nonlinear system (1) exponentially 
converges to that of the stable system (6) if the following dynamic update law 



-ri* + Si 



and the observer system 

x = /(x,u)-*(x,u)[y-ft(x)] 
are adopted. Here a, /3 G R m are as in (20), and fi, £2 £ R m satisfy 



with V m as in (15) and 



and the mapping 



where 



3LV d 2 L 



9u / 9r3u 



r + 



dL 



T 



d 2 L 



3u / dxdu 



<5(x,u) = -e(x,u)P _1 [D/i(x 



*(*,„) > 7±W > o 



cMT 



72 (X) 
with 7i(x, u) > and 72 (x) > defined as in (22) and (23). 



(24) 



(25) 



(26) 



(27) 



(28) 



(29) 



Proof. From the Lyapunov-like function (14), we obtain its time derivative as 



V 



dL\ T d 2 L (dL\ T d 2 L 
d^x) d^ 2 + [dX) dud\ 



dL\ T d 2 L 
— I ^— r 



u + 



dL 



d 2 L 



du I 9A3u 



dL\ T d 2 L . T . 
— I rr^x + e'Pe 



3u/ 3r9u \3u/ 3x3u 
Substituting e in (5), ci and /3 as in (20) and 6 as in (27) into (30), we have 

V = a T ix + p T A + 5 + e T P {/(x,u) - /(x,u) + <t>(x,u)[y - h(x)}} 



(30) 



(31) 



Consider e£S. Since S is convex, according to Mean Value Theorem, there exists mo, mj 6 S 
satisfying 



/(x,u) -/(x,u) = D x f(x + m 1 ,u)e 
y — h(x) = Dh(x + mg)e 



(32) 
(33) 
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Then substituting (24), (26), (32) and (33) into (31), we obtain 

V = -online - p T T 2 p - coV m + e T P [D x /(x + m lr u) + 4>(x,u)D/i(x + m )] e (34) 
After substituting <t>(x, u) as in (28) and R(x, trio) as i n (18)/ (34) can be rewritten as 

V= -a T r 1 a-^ T r 2 ^-a;V,„ + e T PD x /(x + m 1 ,u)e- 6 ^' U ' e T R{x,m )e (35) 

Since the matrices T^ > and T 2 > 0, we have 

V < -ajV m + e T FD x f{x + m lr u)e-^^-e T R{x,m )e (36) 

For e = where x is determined by the observer accurately, we have 

V < -CoV m = -cvV (37) 

Since w > 0, V exponentially converges to zero as t — > +oo. Hence, ( — , 3— I exponentially 

\ a A au J 
converges to zero. 

For any nonzero e G S,letf = r||e|| e. Obviously, v G dS. Then we have 

V < -wV m +i||e|| 2 i/ r PD x /(x + m 1 ,u)i/- ^^||e|| 2 v T R(x,m )i/ (38) 

In the following, we shall show that V converges exponentially to zero for all mg, mj 6 S, 
x G W, u G n, e G S, e ^ and v G dS. 

First let us consider nonzero e G N D S. From v = r || e || e, we have v S M. Since mg, 
mi G S C Q, x G W and u G CI, according to Assumptions 1-5, it follows that 

v T R(x,m )v = (39) 

and 

i/ T PD x /(x + m 1 ,u)i/ < -A: ||i/|| 2 (40) 

with the constant kg > 0. From (38), we have 

V < -wV m - fc ||e|| 2 < -crV (41) 

„ TT fdL dL \ 

with the constant u > 0. Hence, — , — , e I exponentially converges to zero as t — > +oo. 

\ a A ou / 

Then we consider nonzero e G S — N n S, namely, v G dS — M. From (38), taking into account 
(22)-(23), we obtain 

V < -coV m + ^||e|| 2 [ti(x,u) - k r 2 - e{x,u) 72 {x)} (42) 

Since 6(x, u) satisfy the condition (29), we obtain (41) again. Hence, in this case, ( t— -, — , e 

\oA ou 

also exponentially converges to zero as t — > +oo. 
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Since 



dL dL 
dX'du' 



e I exponentially converges to zero as t 



-oo, the closed-loop system 



exponentially converges to 



x = A rf x + B d r 

e = {D x /(x + m 1 ,u)-0(x,u)p- 1 [D/!(x)] T D/i(x + m o )}e 



(43) 



Since A d is a asymptotically stable matrix, we know that x £ Wis bounded. According to 
Assumptions 1 and 4, D x /(x + mj, u), Dh(x) and Dh(x + mrj) are all bounded for mrj, mj £ S 
and u £ Q. From fcrj > 0, we have < 71 (x, u) < +00. According to Assumption 4, we have 
kerR(x,mo) C kerDh(x) which ensures that < i/ T R(x,mrj)v < +00 for every v £ 3S- M, 
rrig 6 S and x 6 W. Thus, we have < 72 (^) < +°°- As a result, < 9(x, u) < +00. From 
(43), we know that e exponentially converges to zero as e exponentially converges to zero. 
Moreover, we have 



A d x - A rf e + B d r 



(44) 



Since e and e exponentially converges to zero, we have the system (1) exponentially converges 
to x = A d x + B d r. This completes the proof. □ 

Consider now the issue of solving (26) with respect to £j and £2- One method to achieve 
a well-defined unique solution to the under-determined algebraic equation is to solve a 
least-square problem subject to (26). This leads to the Lagrangian 



where p £ R is a Lagrange multiplier. The first order optimality conditions 






o, 



dl 
W2 



0, 



'dl 



(45) 



(46) 



leads to the following system of linear equations 



' I m a' 




"ff 




I m P 




li 


= 


_a T /3 T 0_ 




.9 . 









—5 — wV„, 



(47) 



Remark 4. It is noted that Equation (47) always has a unique solution for fj and £2 if any one of ol 
and f> is nonzero. 



4. Example 

Consider the pendulum system 

X2 



X2 

- sin x\ + u\ cos X\ + i(2 sin x\ 
x\+x 2 



(48) 
(49) 
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with x = [x\ X2] £ K. , u = [«i u 2 ] £fl and 

A 



n = |u = [«! « 2 ] T | - 1 < «i < 1, -0.5 < !( 2 < 0.5] 



(50) 



As the system is affine in control and its measurement output y is a linear map of its state x, 
Assumptions 1, 3 and 4 are satisfied automatically. 



Choose 



3 
1 



|ei=-e 2 



For e ^ and e £ ker[l 1], we have e\ = — e 2 and 

e T PD x f{x,u)e\ ei= - e2 

^ x 2 ' — cos*! — U\ sin*! + w 2 cosxi 
= (— cosxi — u\ sinxi + m 2 cos*i + 3)eie 2 | ei =-e 2 

< [— 1.5cos(arctan |) — sin(arctan |) + 3]£i£ 2 | ei= - e2 
= (-1.8028 + 3)e 1 e 2 U=- ez 

= -0.5986||e|| 2 | ei= _ e2 

< - fc ol|e|| 2 | ei =-e 2 

with < kg < 0.5986. Hence, Assumption 5 is satisfied. Let S be the ball of radius r = 1, 
centered at zero and dS is the boundary of S. Define M C 3S and 

M= lv= [vj v 2 ] T G K 2 : |M| = 1, 3vjv 2 + 1.8028|i/iv 2 | < -fc } 
Obviously, 

dS - M = {1/ = [vj i/ 2 ] T G R 2 : ||v|| = 1, 3i/!i/ 2 + 1.8028|v 1 v 2 | > -k Q \ 
As 7! (x, «) = 3 x 1.8028 + fc and 

7 2 (x) =min|(i/ 1 +i/ 2 ) 2 , vg3S-m]=1 
?i(x,u) 



choosing A:g = 0.5, we have 



72(x) 



3 - 1.8028 
35.8699. Let 6(x,u) = 36 > 35.8699 and we have 



4>(x,u) = -[12 36] T . 

Now the nonlinear observer becomes 

*2 



-1'2 



— sin *i + Mj cos X\ + m 2 sin Xj 
Choose the reference model (6) where 



+ 



(y-x 1 -x 2 ) 



1 
-25 -10 



B, 
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and the reference is given by 



15 



10 



t- t 2 
tf-tz 



15 



t f -t 2 



+ 10 



t- 1 2 
tf-h 



o<t<t 1 
t x <t<t 2 

+ rp t 2 <t<t f 
t > t f 



with t\ = 10s, t 2 = 20s, tf = 30s and re = 0.5. Obviously, Assumption 2 is satisfied. 

Set Hx = 0, H 2 = 10~ 4 I 2 , a; = l,Ti = T 2 = 2I 2 , and xi(0) = 0.3 and x 2 (0) = 0.5. Using the 
proposed approach, we have the simulation result of the pendulum system (48)-(50) shown in 
Figures 2-5 where the control u 2 is stuck at —0.5 from t = 12s onward. 

From Figure 2, it is observed that the estimated states X\ and x 2 converge to the actual 
states X\ and x 2 and match the desired states X\£ and Xyj 'well, respectively, even when 
u 2 is stuck at —0.5. This observation is further verified by Figure 3 where both the state 
estimation errors e\{= X\ — X\) and e 2 (= x 2 — x 2 ) of the nonlinear observer as in (4) and 
the matching errors Tj(= 0) and t 2 (= — sinfj + U\ cosfi + «2 s i n ^l + 25fi + 10^2 — 25r) as 
in (8) exponentially converge to zero. Moreover, Figure 4 shows that the control U\ roughly 
satisfies the control constraint U\ £ [—1,1] while the control u 2 strictly satisfies the control 
constraint u 2 6 [—0.5,0.5]. This is because, in this example, the Lagrange multiplier Aj is first 
activated by the control U\ < — 1 at t = (see Figure 5 where Aj is no longer zero from t = 0), 
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Fig. 2. Responses of the desired, estimated and actual states 
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Fig. 3. Responses of estimation error and matching error 
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Fig. 4. Responses of control u 

and then the proposed dynamic update law forces the control U\ to satisfy the constraint 
U\ £ [—1, 1]. It is also noted from Figure 5 that the Lagrange multiplier A 2 is not activated in 
this example as the control 112 is never beyond the range [—0.5, 0.5] . In addition, the output y 
and the Lyapunov-like function V m are shown in Figure 6. From Figure 6, it is observed that 
the Lyapunov-like function V m exponentially converges to zero. 
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Fig. 5. Responses of Lagrangian multiplier A 
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Fig. 6. Responses of output y and Lyapunov-like function V m 

5. Conclusions 

Sufficient Lyapunov-like conditions have been proposed for the control allocation design via 
output feedback. The proposed approach is applicable to a wide class of nonlinear systems. 
As the initial estimation error e(0) need be near zero and the predefined dynamics of the 
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closed-loop is described by a linear stable reference model, the proposed approach will 
present a local nature. 
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1. Introduction 

Flexible-link robotic manipulators have many advantages with respect to conventional rigid 
robots. They are built by lighter, cheaper materials, which improve the payload to arm 
weight ratio, thus resulting in an increase of the speed with lower energy consumption. 
Moreover, due to the reduced inertia and compliant structure, these lightweight arms can be 
operated more safely and are more applicable for the delicate assembly tasks and interaction 
with fragile objects, including human beings. 

The control for robot manipulators is to determine the time history of joint inputs to cause 
the end-effector to execute a commanded motion. There are many control techniques and 
methodologies that can be applied to the control of the manipulators. The specific control 
method and its implementation ways can have a significant impact on the performance of 
the manipulator and consequently on the range of its possible applications. In addition, the 
mechanical design of the manipulator itself will influence the type of control scheme 
needed. However, in order to improve the control performance, more sophisticated 
approaches should be found. 

The control for flexible joint system has attracted a considerable amount of attention during 
the past few years. There are PD, inverse dynamics, and the force control approach for the 
feedback control strategies of flexible joint manipulator. (1989, MARK W. SPONG), an 
integral manifold approach to the feedback control of flexible joint robots (1987, MARK W. 
SPONG, KHASHAYAR KHORASANI, and PETAR V. KOKOTOVIC), and the nonlinear 
feedback control of flexible joint manipulators: a single link case study, (1990, K. 
KHORASANI). The basic idea of feedback linearization is to construct a nonlinear control 
law as a so-called inner loop control which, in the ideal case, exactly linearizes the nonlinear 
system after a suitable state space change of coordinates. The designer can design a second 
stage or outer loop control in the new coordinates to satisfy the traditional control design 
specifications such as tracking, disturbance rejection, and so forth. Since the feedback 
linearization of flexible joint manipulator is a fourth order integrator system, so we 
proposed a three stage design method, the first is nonlinear feedback to get integrator 
system, the second is pole placement to get expect performance, and the third is to use PFC 
to reject disturbance and uncertainty, since they can not be exactly cancelled by nonlinear 
feedback, coupling effects of the joint flexibility. More accurate description of robot 
dynamics may include fast actuator dynamics and joint-link flexibility, and so on. 
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2. Equations of motion 

Consider the single-link arm shown in Figure 1 consisting of a flexible joint. 




Fig. 1. Single-link robot with joint flexibility 

The kinetic energy of the manipulator is a quadratic function of the vector q 



1 lA 



'/; 



(1) 



where the nxn inertia matrix D(q) is symmetric and positive definite for each q e s Jt" . 



The potential energy V = V(q) is independent of ij . We have remarked that robotic 
manipulator satisfies this condition. 



The Euler-Lagrange equations for such a system can be derived as follows. Since 

L = K-V = \f j d lj (q)q l q j -V(q) 



(2) 



we have 



and 



Also 



% j 



d dL 
dt dq k 



d_ 
'dt 



,ddi. 
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% 2 hj dq k ' ' dq k 
Thus the Euler-Lagrange equations can be written as 

^ —,{8d ]d iddy] dV 

E^)^E{^-^4 i ---,^V-,« (3) 

By interchanging the order of summation and taking advantage of symmetry, we can show 
that 

dd kj\.. 1 v I dA n dd ki 



Hence 



M\n=m^\n 



Tj { fyi 2 % J tj 2 { d li d ij 5( ?/t 



The term 



1 8d kj dd u dd tj 

are known as Christoffel symbols. Note that, for a fixed k, we have c,-. t = c ft , which reduces 
the effort involved in computing these symbols by a factor of about one half. Finally, if we 
define 

dV «x 

Pk = t— (5) 

then the Euler-Lagrange equations can be written as 

Z^(^/ + Z C »;*(?^; + %W = T *' fc = V",n (6) 

In the above equation, there are three types of terms. The first involve the second derivative 
of the generalized coordinates. The second are quadratic terms in the first derivatives of q, 
where the coefficients may depend on q. These are further classified into two types. Terms 
involving a product of the type q i are called centrifugal, while those involving a product of 
the type <j;4; where i^j are called Coriolis terms. The third type of terms are those 
involving only q but not its derivatives. Clearly the latter arise from differentiating the 
potential energy. It is common to write (6) in matrix form as 

D(q)q + C(q,q)q + g(q) = T (7) 
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where the k, j-th element of the matrix C(q,q) is defined as 

c kj =Zc p (<?)4,. = £- J_i + -Ji-.--l 
3. Feedback linearization design for inner loop 

We first derive a model similar to (6) to represent the dynamics of a single link robot with joint 
flexibility. For simplicity, ignoring damping of the equations of motion, system is given by 

D (<7i)3i + M<h'<7i)<7i +K(<h -<? 2 ) = ( 8 ) 

Jq 2 -K(q 1 -q 2 ) = u (9) 

In state space, which is now 9? " , we define state variables in block form 



1 — W1 Xi — Cf-i 

X 3 = 12 X 4 = tfl 



(10) 



Then from (8)- (9) we have 

x 1= x 2 (U) 

x 2 =-D(x 1 )- 1 {h(x 1 ,x 2 ) + K(x 1 -x 3 )} (12) 

i 3 =x 4 (13) 

x 4 = r 1 K(x 1 -x 3 ) + r 1 M (14) 

This system is then of the form 

x = f(x) + G(x)u (15) 

For a single-input nonlinear system, /(x) and ^(x) are smooth vector fields on s Jt" , 
/(0) = , and iieSR, is said to be feedback linearizable if there exists a region U in 
9i" containing the origin, a diffeomorphism T: L7 — » 9?" , and nonlinear feedback 

u = a(x) + p(x)v (16) 

with jB(x) * on U, such that the transformed variables 

y = T(x) (17) 

satisfy the system of equations 

if = Ay + bv (18) 

where 
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10 
1 







1 





In the single-link case we see that the appropriate state variables with which to define the 
system so that it could be linearized by nonlinear feedback on the link position, velocity, 
acceleration, and jerk. Following the single-input case, then, we can apply same action 
on the multi-link case and derive a feedback linearizing transformation blockwise as 
follows, 



y 1 =T 1 (x) = X 1 

y 2 =T 2 (x) = y 1 =x 2 

1/3 = T 3 {x) = y 2 =x 2 = -D{XiY l {h(x 1 ,x 2 ) + K(* a -x 3 )} 
y 4 = T 4 (x) = y 3 = -4:[D(x l r 1 ]{h(x 1 ,x 2 ) + K(x, -x 3 )}- Dfo)" 1 ^* 



dt 



ex 



cli 



(19) 
(20) 

(21) 
(22) 



+— l-D(Xi) 1 (h(x 1 ,x 2 ) + K(x t - x 3 ))] + K(x 2 - x 4 )} = a 4 (x 1 ,x 2 ,x 3 ) + D(x x ) a Kx 4 

dx 2 

where for simplicity we define the function a 4 to be that in the definition of y 4 except the 
last term, which is D~ Kx 4 . Note that x 4 appears only in this last term so that a 4 depends 
only on x 1 ,x 2 ,x 3 . 

As in the single-link case, the above mapping is a global diffeomorphism. Its inverse can be 
found by 



X 1 =y 1 

x 2 =y 2 

x 3 =y 1 +K- 1 {D{y 1 )y 3 +h{y 1 ,y 2 )) 

x 4 = K- 1 D(y 1 )(y i - a i {y 1 ,y 2 ,y 3 )) 
The linearizing control law can now be found from the condition 

y A = v 



(23) 
(24) 
(25) 
(26) 

(27) 



where v is a new control input. Computing y 4 from (22) and suppressing function 
arguments for brevity yields 
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ca. 



da. 



dx x dx 2 

■ a(x) + b(x)u 



D-\h + K( Xl - x 3 )) + '^-x i + — [D- 1 ]Kx i + D- 1 K{r 1 K(x 1 -x 3 ) + }' r u) 



ox 



dt 



(28) 



where 



a(x)-- 



uu 



x 2 --^D- 1 (h + K(x 1 -x 3 )) 



da a 



dx 1 dx 2 ' dx 3 dt 



+ ^[D- 1 ]Kx 4 +D- 1 KJ- 1 K(x 1 -x 3 ) (29) 



b(x) = D~ 1 (x)K}~ 1 u 
Solving the above expression for u yields 

u = b(x)~ 1 (v-a(x)) 

=: a(x) + P(x)v 

where p{x) = ]K^D(x) and a{x) = -b{x)^ a(x) 



(30) 

(31) 
(32) 



10 0" 




"0" 


10 







1 


y + 










I 



With the nonlinear change of coordinates (19)-(22) and nonlinear feedback (32) the 
transformed system has the linear block form 



(33) 



=: Ay + bv 

where I = nxn identity matrix, = nxn zero matrix, y = (y 1 ,y 2 ,y 3 ,J/ 4 ) e 9? ", and 
»e*R". The system (33) represents a set of n decoupled quadruple integrators. 

4. Outer loop design based on predictive function control 

4.1 why use predictive function control 

The technique of feedback linearization is important due to it leads to a control design 
methodology for nonlinear systems. In the context of control theory, however, one should 
be highly suspicious of techniques that rely on exact mathematical cancellation of terms, 
linear or nonlinear, from the equations defining the system. 

In this section, we investigate the effect of parameter uncertainty, computational error, 
model simplification, and etc. We show that the most important property of feedback 
linearizable systems is not necessarily that the nonlinearities can be exactly cancelled by 
nonlinear feedback, but rather that, once an appropriate coordinate system is found in 
which the system can be linearized, the nonlinearities are in the range space of the input. 
This property is highly significant and is exploited by the predictive function control 
techniques to guarantee performance in the realistic case that the nonlinearities in the 
system are not known exactly. 
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Consider a single-input feedback linearizable system. After the appropriate coordinate 
transformation, the system can be written in the ideal case as 

1/1=1/2 

■ (34) 

y„= v = fi~ 1 (x)[u-a(x)] 

provided that u is given by (16) in order to cancel the nonlinear terms a(x) and j6(x) . 

In practice such exact cancellation is not achievable and it is more realistic to suppose that 
the control law u in (16) is of the form 

u = a(x) + j3(x)v (35) 

where d(x) , f)(x) represent the computed versions of a(x) , /3(x) , respectively. These 
functions may differ from the true a(x) , /3(x)for several reasons. Because the inner loop 
control u is implemented digitally, there will be an error due to computational round-off 
and delay. Also, since the terms a(x) , /3(x) are functions of the system parameters such as 
masses, and moments of inertia, any uncertainty in knowledge of these parameters will be 
in reflected in a(x) , P{x) . In addition, one may choose intentionally to simplify the control 
u by dropping various terms in the equations in order to facilitate on-line computation. If 
we now substitute the control law (35) into (34) we obtain 

1/1 = 1/2 

: (36) 

tin = V = /^ (*)[«(*) + P{ x ) v - «(*)] 

= v + 7j(y 1 ,---,y„,v) 

where the uncertainty 77 is given as 

l{yx,-,y n ,v) = {{p- l p-l)v + p-\a-a))\ y= ^ {x) (37) 

The system (36) can be written in matrix form as 

y = Ay + b{v + r 1 {y,v)} (38) 

where A and b are given by (18). For multi-input case, similar to (33), and if v e iR m , and 

rj : 9?" x 9? m — > <R m . Note that the system (38) is still nonlinear whenever 77 ^0 . The practical 
implication of this is solved by the outer loop predictive function control (PFC). 

The system (38) can be represented by the block diagram of Figure 2. The application of the 
nonlinear inner loop control law results in a system which is "approximately linear". A 
common approach is to decompose the control input v in (38) into two parts, the first to 
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stabilize the 'nominal linear system' represented by (38) with 77 = . In this case v can be 
taken as a linear state feedback control law designed to stabilize the nominal system and/ or 
for tracking a desired trajectory. A second stage control Av is then designed for robustness, 
that is, to guarantee the performance of the nominal design in the case that r] ^ . Thus the 
form of the control law is 



u = d{x) + J3{x)v 
v = -Ky + Av 
Av = PFC(y r ) 



(39) 
(40) 
(41) 



where Ky is a linear feedback designed to place the eigenvalues of A in a desired location, 
Av represents an additional feedback loop to maintain the nominal performance despite the 
presence of the nonlinear term tj . y r is a reference input, which can be chosen as a signal for 
tracking a desired trajectory. 
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Fig. 2. block diagram for PFC outer loop design 



4.2 Predictive function control 

All MPC strategies use the same basic approach i.e., prediction of the future plant outputs, 
and calculation of the manipulated variable for an optimal control. Most MPC strategies are 
based on the following principles: 

Use of an internal model 

Its formulation is not restricted to a particular form, and the internal model can be linear, 
nonlinear, state space form, transfer function form, first principles, black-box etc. In PFC, 
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only independent models where the model output is computed only with the present and 
past inputs of the process models are used. 

Specification of a reference trajectory 

Usually an exponential. 

Determination of the control law 

The control law is derived from the minimization of the error between the predicted output 
and the reference with the projection of the Manipulated Variable (MV) on a basis of functions. 

Although based on these principles the PFC algorithm may be of several levels of 
complexity depending on the order and form of the internal model, the order of the basis 
function used to decompose the MV and the reference trajectory used. 

4.3 First order PFC 

Although it is unrealistic to represent industrial systems by a first order system, as most of 
them are in a higher order, some well behaved ones may be estimated by a first order. The 
estimation will not be perfect at each sample time, however, the robustness of the PFC will 
help to maintain a decent control. 

If the system can be modelled by a first order plus pure time delay system, then the 
following steps in the development of the control law are taken. 

Model formulation 

In order to implement a basic first order PFC, a typical first order transfer function equation 
(42) is used. 

y M (s) = -^-u(s) (42) 

T M S + 1 

Note that the time delay is not considered in the internal model formulation and in this case 
K M is equal to one. The discrete time formulation of the model zero-order hold equivalent 
is then obtained in (43). 

y M (k) = ay M (k-l) + K M (l-a)u(k-l) (43) 

T 
where a = exp( — ) . If the manipulated variable is structured as a step basis function: 

I'm 

y L (k + H) = a H y M (k) (44) 

y F (k + H) = K M (l-a H )u(k) (45) 

Where, y L and y F are respectively, the free (autoregressive) and the forced response of y M . 
Reference trajectory formulation 
If y R is the expression of the reference trajectory, then at the coincidence point H: 
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C(k + H)-y R (k + H) = A H (C(k)-y P (k)) 



(46) 



thus: 



y R (k + H) = C(k)-A H (C(k)-y P (k)) 



(47) 



Predicted process output 



The predicted process output is given by the model response, plus a term given the error 
between the same model output and the process output: 



y p {k + H) = y M (k + H) + (y p (k) - y M {k)) 

where y M (k + H) = y L (k + H) + y F (k + H) = a H y M (k) + K M (l-a H )u(k) . 
Computation of the control law 
At the coincidence point H: 

y R (k + H) = y P (k + H) 
Combining (44), (45), (47) and (48) yields 

C(k)-A H (C(k)-y P (k))-y P (k) = y M (k + H)-y M (k) 
Replacing y M (k + H) by its equivalent in equations (44) and (45) we obtain: 

C(k)(l-A H )-y P (k)(l-A H ) + y M (k)(l-a H ) = K M (l-a H )u(k) 
Solving for u(k) the final result is the control law given in (52). 

{C{k)-y P {k)){l-A R )y M {k) 



u{k) 



(48) 



K M (l-a H ) 



(49) 



(50) 



(51) 



(52) 



4.4 Case of a process with a pure time delay 

In the linear case, a process with a pure time delay can be expressed in terms of a delay -free 
part, plus a delay added at the output, as in Fig. 3. 



yPdelay 




Fig. 3. Process with time delay 

The value y P & e i m at time k is measured, but not y p . In order to take into account the delay 
in a control law formulation, prior knowledge of the delay value d is needed. y p can be 
estimated as: 
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y P ( k ) = ypddayi k ) + y M ( k ) -Vm^- d ) ( 53 ) 

4.5 Tuning in PFC 

According to the three principles of PFC, tuning is a function of the order of the basis 
constructing the MV, the reference trajectory, the control horizon and the CLRT value. 

The influence of the PFC parameters is given in Table 1, where the influence of various 
PFC parameters is on precision (Steady State Resp.), transient response and robustness 
are graded between (indicating minimum influence) and 1 (indicating maximum 
influence). 

SS Resp. Transient Resp. Robustness 

Basis function 10 

Reference trajectory 1 1/2 

Coincidence horizon 1/2 1 

Table 1. Effect of PFC parameters in tuning 

In most cases, an exponential reference trajectory is chosen along with a single coincidence 
horizon point (H = 1) and a zero order basis function (Richalet, 1993). Considering the 
known Open Loop Response Time of the system (OLRT), one can choose the CLRT value 
given by the ratio OLRT/ CLRT. This ratio then becomes the major tuning parameter 
shaping the system output and MV, dictating how much overshoot occurs and ensuring 
stability, on the condition that the internal model is accurate enough. For slow processes, 
e.g., heat exchange systems, a ratio of 4 or 5 is found most suitable, and ensures a stable 
PFC. 

5. Simulation 

Consider the single link manipulator with flexible joint shown in Figure 1. Choosing q 1 and 
q 2 as generalized coordinates, the kinetic energy is 

K = 1^+1^2 (54) 



The potential energy is 



V = MgL(l-cosq 1 ) + h(q 1 -q 2 ) 2 (55) 



The Lagrangian is 



L = K-V=hq 2 1+ ±Jq 2 2 -MgL(l-cosq 1 )-h(q 1 -q 2 ) 2 (56) 



Therefore we compute 
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dL T . dL T . 

T7- = 111 T7- = 7<?2 






-— = -MgLsin((j 1 )-% 1 -(j 2 ) — = %i-?2) 

Therefore the equations of motion, ignoring damping, are given by 

Iq\ + MgL sinfijj ) + k(q- i - q 2 ) = 

JiJ2-fc(<7i-<7 2 ) = M 



(57) 
(58) 
(59) 

(60) 
(61) 



Note that since the nonlinearity enters into the first equation the control u cannot simply be 
chosen to cancel it as in the case of the rigid manipulator equations. In other words, there is 
no obvious analogue of the inverse dynamics control for the system in this form. 



In state space we set 



1 — 1 1 2 — t1 



x 3 — q 2 x A -q 2 



and write the system (60)- (61) as 



MgL . k. 

x 2 = — sm(x 1 )-y(x 1 -x 3 ) 



x 4=y( x i- x 3) + y« 



The system is thus of the form (15) with 



/(*) = 



MgL . k. 

■—f-sm(x 1 )--(x 1 -x 3 ) 



J 



{x 1 -x 3 ) 



g(x)- 



(62) 



(63) 



Therefore n=4 and the necessary and sufficient conditions for feedback linearization of this 
system are that 
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rankig,ad f (g),ad 2 f (g),ad 3 f (g)\ = rank 
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1 




k 






u 




u 



(64) 



/ / 2 

which has rank 4 for k > , I,/ < co . Also, since vector fields lg,adf(g),ad r{g)\ are constant, 
they form an involutive set. 



{#/«<*/ (sW/(g)} 



(65) 



To see this it suffices to note that the Lie Bracket of two constant vector fields is zero. Hence 
the Lie Bracket of any two members of the set of vector fields in (65) is zero which is trivially 
a linear combination of the vector fields themselves. It follows that the system (60)- (61) is 
feedback linearizable. The new coordinates 



y,=T, i = i,-A 

are found from the conditions (67)- (68) 

<dT lr ad k f (g)>=0 fc = 0,l,---,n-2 



<dT 1 ,ad n f 1 (g)>*0 



(66) 

(67) 
(68) 



with n=4, that is 



<dT lr g >=0 
<dT 1 ,[f,g]>=0 
<dT lr adl(g) >=0 

<dT lt ad}(g)>*0 
Carrying out the above calculations leads to the system of equations 



OX j 



= 



dT ^-o ^ = 



dx. 



dx. 



and 



dx i 



-*0 



(70) 
(71) 
(72) 

(73) 
(74) 
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From this we see that the function T a should be a function of x 1 alone. Therefore, we take 
the simplest solution 



and compute 



y 1 =T 1 =x 1 



y 2 = T 2 =< dT-y,f >= x 2 



y 3 = T 3 =< dT 2 ,f >= ^-sin(x 1 ) -j( Xl -x 3 ) 



y 4 =T 4 =<dT 3 ,f>= ^-cos(x 1 )xx 2 -y(x 2 -x 4 ) 

The feedback linearizing control input u is found from the condition 

1 



(75) 

(76) 
(77) 

(78) 



<dT i/g > 



(v-<dT A ,f >) 



II 



(79) 



(v - fl(x)) = fS(x)v + a(x) 



where 



i \ M S L ■ i w 2 MgL k k k k MgL 

« (*) : = — j-sm( Xl )(x 2 + — ^cos(x 1 ) + y) + y (xj - x 3 )(y + y + — ^-< 



s(*i)) (80) 



Therefore in the coordinates y 1 ,--,i/ 4 with the control law (79) the system becomes 

ill =Vi J/2 = Vi 

y 3 =yi iii = v 



or, in matrix form, 



where 



if = Ay + bv 



(81) 



10 0" 




"0" 


10 
1 


b = 











1 



The transformed variables J/i,---,i/4 are themselves physically meaningful. We see that 

y 1 = x 1 =link position 
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y 2 = x 2 =link velocity 
3/3 = Vi = link acceleration 

y 4 =y 3 =link jerk 

Since the motion trajectory of the link is typically specified in terms of these quantities they 
are natural variables to use for feedback. 

For given a linear system in state space form, such as (81), a state feedback control law is an 
input v of the form 



v = -k T y + r = -£%,•• 

1=1 



(82) 



where fc,- are constants and r is a reference input. If we substitute the control law (82) into 
(81), we obtain 



y = (A-bk )y + br 



(83) 



Thus we see that the linear feedback control has the effect of changing the poles of the 
system from those determined by A to those determined by A-bk 

When the parameters are chosen fc x = 62.5 , k 2 = 213.6 , k 3 = 204.2 A: 4 = 54 , we can get step 
responses in Figure 4. where kl, k2, k3 and k4 are linear feedback coefficients to place the 
eigenvalues of A in a desired location. 






1 











1 











62.5 


-213.8 


-204.2 



" 


Vi 




Tf 





3/2 


+ 





1 


1/3 







-54 


1/4 




1 



(84) 



y = [l 0] 



(85) 



The internal model parameter: K M = 0.016 , T M = 3 , d = 8 , and the coincidence point H=10. 
Response time of reference trajectory is 0.01, and sample time is 0.01. 
For the uncertainty 77 , the system (38) can be written in matrix form as 

y = {A-bk T )y + b{r + i 1 {y,v)} 

then use predictive function control strategy to reduce or overcome uncertainty of nonlinear 
feedback error tj(y,v) ,and simulation result is shown in Figure 5 for rj(y,v) = 10%y r . 
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Fig. 4. link position output 
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Fig. 5. link position output with uncertainty rejection 
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6. Conclusion 

A new three stage design method is presented for the single link manipulator with flexible 
joint. The first is feedback linearization; the second is to use pole placement to satisfy 
performance, and the third is to develop predictive function control to compensate 
uncertainty. Finally, for the same uncertainty, robustness is better than traditional method. 
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1. Introduction 

Several classes of general hybrid and switched dynamic systems have been extensively 
studied, both in theory and practice [3,4,7,11,14,17,19,26,27,30]. In particular, driven by 
engineering requirements, there has been increasing interest in optimal design for hybrid 
control systems [3,4,7,8,13,17,23,26,27]. In this paper, we investigate some specific types of 
hybrid systems, namely hybrid systems of mechanical nature, and study the corresponding 
hybrid OCPs. The class of dynamic models to be discussed in this work concerns hybrid 
systems where discrete transitions are being triggered by the continuous dynamics. The 
control objective (control design) is to minimize a cost functional, where the control 
parameters are the conventional control inputs. 

Recently, there has been considerable effort to develop theoretical and computational 
frameworks for complex control problems. Of particular importance is the ability to operate 
such systems in an optimal manner. In many real-world applications a controlled mechanical 
system presents the main modeling framework and can be specified as a strongly nonlinear 
dynamic system of high order [9,10,22]. Moreover, the majority of applied OCPs governed 
by sophisticated real-world mechanical systems are optimization problems of the hybrid 
nature. The most real-world mechanical control problems are becoming too complex to allow 
analytical solution. Thus, computational algorithms are inevitable in solving these problems. 
There is a number of results scattered in the literature on numerical methods for optimal 
control problems. One can find a fairly complete review in [3,4,8,24,25,29]. The main idea 
of our investigations is to use the variational structure of the solution to the specific two point 
boundary-value problem for the controllable hybrid-type mechanical systems in the form of 
Euler-Lagrange or Hamilton equation and to propose a new computational algorithm for the 
associated OCR We consider an OCP in mechanics in a general setting and reduce the initial 
problem to a constrained multiobjective programming. This auxiliary optimization approach 
provides a basis for a possible numerical treatment of the original problem. 

The outline of our paper is as follows. Section 2 contains some necessary basic facts 
related to the conventional and hybrid mechanical models. In Section 3 we formulate 
and study our main optimization problem for hybrid mechanical systems. Section 4 deals 
with the variational analysis of the OCP under consideration. We also briefly discuss the 
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computational aspect of the proposed approach. In Section 5 we study a numerical example 
that constitutes an implementable hybrid mechanical system. Section 6 summarizes our 
contribution. 

2. Preliminaries and some basic facts 

Let us consider the following variational problem for a hybrid mechanical system that is 
characterized by a family of Lagrange functions {Lp ; }, p,- 6 V 



Joti' 1 '- 1 "" (It 

subject to q{0) = Cq, q(\) = C\, 



f 
minimize / £ %_ lA )(0 L P( ( f ' "7(0/ 4(0)* 



where T 5 is a finite set of indices (locations) and q(-) (q(t) G R") is a continuously differentiable 
function. Here pu., f ,\ (■) are characteristic functions of the time intervals [f;_i, f;), z = 1, ..., r 
associated with locations. Note that a full time interval [0, 1] is assumed to be separated into 
disjunct sub-intervals of the above type for a sequence of switching times: 

r:={t = 0,h,...,t r = l}. 

We refer to [3,4,7,8,13,17,23,26,27] for some concrete examples of hybrid systems with the 
above dynamic structure. Consider a class of hybrid mechanical systems that can be 
represented by n generalized configuration coordinates q\,...,q n . The components q\{t),A = 
1, ..., n of q(t) are the so-called generalized velocities. Moreover, we assume that L pi (t, ■,•) 
are twice continuously differentiable convex functions. It is well known that the formal 
necessary optimality conditions for the given variational problem (1) describe the dynamics of 
the mechanical system under consideration. This description can be given for every particular 
location and finally, for the complete hybrid system. In this contribution, we study the hybrid 
dynamic models that free from the possible external influences (uncertainties) or forces. The 
optimality conditions for mentioned above can be rewritten in the form of the second-order 
Euler-Lagrange equations (see [1]) 

A dL pi (t,q,q) dl pi (t,q,q) 

dt dq A dq A ' '•"' ' (2) 

q{0) = c , q(l) = c\, 

for all pi G V '. The celebrated Hamilton Principle (see e.g., [1]) gives an equivalent variational 
characterization of the solution to the two-point boundary -value problem (2). 

For the controllable hybrid mechanical systems with the parametrized (control inputs) 
Lagrangians L„ j (t, q, q, u), j)j G V we also can introduce the corresponding equations of 
motion 

d dL Vi (t,q,q,u) dL pi (t,q,q,u) 

dt dq A dq A ~ ' (3) 

q(Q) = Cq, q(l) = c lr 
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where «(•) £ U is a control function from the set of admissible controls U. Let 

U := {u £ R m : \ v < u v < b 2/V , v = l,-,m}, 
W:=K)eL 2 ffl ([0,l]) : v(t) e LI a.e. on [0,1]}, 

where b\ v , b^.v, v = l,—,tn are constants. The introduced set U provides a standard example 
of an admissible control set. In this specific case we deal with the following set of admissible 
controls U nC„(0, 1). Note that L Pi depends directly on the control function «(•). Let us 
assume that functions L pi (t, •, -,u) are twice continuously differentiable functions and every 
L pi (t, q, q, •) is a continuously differentiable function. For a fixed admissible control «(■) we 
obtain for all p,- £ V the above hybrid mechanical system with L p; (f, g, 4) = £-p ! (^<?/4' M W)- 
It is also assumed that L„.(t, q,-, u ) are strongly convex functions, i.e., for any 



(t,q,q,u) G R x R" x R" x R m , £ e R" 



the following inequality 



" d 2 L Pi (t,q,q,u) " 3 

holds for all p, G P. This natural convexity condition is a direct consequence of the 
classical representation for the kinetic energy of a conventional mechanical system. Under 
the above-mentioned assumptions, the two-point boundary-value problem (3) has a solution 
for every admissible control u(-) £ U [18]. We assume that (3) has a unique solution for every 
m(-) £ U. For an admissible control w() £ U the solution to the boundary -value problem (3) 
is denoted by q u {-)- We call (3) the hybrid Euler-Lagrange control system. Let us now give a 
simple example of the above mechanical model. 

Example 1. We consider a variable linear mass-spring system attached to a moving frame that is a 
generalization of the corresponding system from [22]. The considered control u(-) £ U f|C{(0, 1) is 
the velocity of the frame. By a>p t we denote the variable masses of the system. The kinetic energy 

K = 0.5co pi (q + u) 2 

depends on the control input «(■). Therefore, we have 

L pi (q,q,u) = 0.5(cOp i (q + u) 2 — nq 2 ), k £ R + 

and 

d dh Vi {t,q,q,u) dh pi {t,q,q,u) ,..,.,, n 

It Iq dq =<"*(* + «0+*f = 0. 

By k we denote here the elasticity coefficient of the spring system. 

Note that some important controlled mechanical systems have a Lagrangian function of the 
following form (see e.g., [22]) 

m 

L pi (t, q,q,u) = L° p .(t, q,q)+ ^ q v u v . 

v=\ 
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In this special case we easily obtain 

d dL°.(t,q,q) dL°.(t,q,q) ( u A A = l,...,m, 



dt dqx dq\ [0 \ = m + l,...,n. 

Note that the control function u() is interpreted here as an external force. 

Let us now consider the Hamiltonian reformulation for the controllable Euler-Lagrange 
system (3). For every location p; from V we introduce the generalized momenta 

s A := L pi {t,q,q,n)/dq A 

and define the Hamiltonian function H„ t (t, q, s, u) as a Legendre transform applied to every 
L„ ( (t,q r q,u), i.e. 

n 
H pj {t,q,s,u) := £ s A q A - L pi (t,q,q,u). 
A=l 

In the case of hyperregular Lagrangians L„ j (t, q, q, u) (see e.g., [1]) the Legendre transform, 
namely, the mapping 

Cp t : (t,q,q,u) — > (t,q,s,u), 

is a diffeomorphism for every p,- E V ,. Using the introduced Hamiltonian H(t, q, s, u), we can 
rewrite system (3) in the following Hamilton-type form 

, % dH„.(t,q,s,u) 

u <t \ (4) 

Under the above-mentioned assumptions, the boundary-value problem (4) has a solution for 
every «(•) G W. We will call (4) a Hamilton control system. The main advantage of (4) in 
comparison with (3) is that (4) immediately constitutes a control system in standard state 
space form with state variables (q,s) (in physics usually called the phase variables). Consider 
the system of Example 1 with 

H pi (q,s,u) = ^co p; (q 2 - u 2 ) + -Kq 2 - su. 

The Hamilton equations in this case are given as follows 

dH Pl (q,s,u) 1 

q = ! -^r = s — u, 

ds w pi 

dH p (q,s,u) 

i = Tq = -' Cq - 

Clearly, for 

m 
L p , (£, q, q, u) = L° Pi (t,q,q)+ J^ q v , u v 

v=\ 
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we obtain the associated Hamilton functions in the form 

m 

H pi (t,q,s,u) = H°.(t,q,s) - J^ q v u v , 

v=\ 

where h2. (t,q,s) is the Legendre transform of L p At, q,q). 

3. Optimization of control processes in hybrid mechanical systems 

Let us formally introduce the class of OCPs investigated in this paper: 

minimize / := jf £ ^ ti _ lA) {t)f°.(q u (t),u(t))dt 
subject to u(t) elite [0,1], t t 6 x, i= l,...,r, 

where fS. : [0, 1] x R" X R m — > R be continuous and convex on R" x R m objective functions. 
We have assumed that the boundary-value problems (3) have a unique solution q % (•) and that 
the optimization problem (5) also has a solution. Let (q°P (•), M°P f (•)) be an optimal solution 
of (5). Note that we can also use the associated Hamiltonian-type representation of the initial 
OCP (5). We mainly focus our attention on the application of direct numerical algorithms to 
the hybrid optimization problem (5). A great amount of works is devoted to the direct or 
indirect numerical methods for conventional and hybrid OC problems (see e.g., [8,24,25,29] 
and references therein). Evidently, an OC problem involving ordinary differential equations 
can be formulated in various ways as an optimization problem in a suitable function space 
and solved by some standard numerical algorithms (e.g., by applying a first-order methods 
[3,24]). 

Example 2. Using the Euler-Lagrange control system from Example 1, we now examine the 
following OCP 

a >■ 

minimize/ :=- / Jj ^[t,-_i,t,-) ( f )* : p<- ("(*) + "?( f )) df 
Jo l=1 

subject to q(t) -\ q(t) = —u(t), i = 1, ..., r, 

q (0) = Q, q(l) = 1, 

«(•) ec}(0,i), o< u(t) <ivf e [o,i], 

where kp t are given (variable) coefficients. Let Wp j > 4k /n 2 for every pi e V. The solution (?"(•) of 
the associated boundary-value problem can be written as follows 

q u {t) = Cf sin (tJic/ujp,) - / Jk/co P{ sin(JK/cj pi (t - l))u{l)dl, 

where t e [f,_i, f,), i = 1, ..., r and 

Cf= J- [1+ / J K /co Pi sm{VU^{t - l))u(l)dl] 

sm^K/cOp i Jo V r 
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is a constant in every location. Consequently, we have 

J = - [t %_ lA ) (0**, [«(0 + ?*(*)]# = 

jo i=1 

- I LPlt^,t,)(t)kpMt)+C?sm(tJK/cv pi ) 
jo i=1 

— / Jk/ tOp i sin (J K/Wp.(t — l))ii(l)dl]dt. 

Let now k Pi = 1 for all p, £ V. Using the Hybrid Maximum Principle (see [4]), we conclude that 
the admissible control «°P'(£) = 0.5 is an optimal solution of the given OCP. Note that this result is 
also consistent with the Bauer Maximum Principle (see e.g., [2] ). For m°P'(-) we can compute the 
corresponding optimal trajectory given as follows 

r'(0 = ™ ( '^ ,*e ft-!, *,•),.■ = ! r. 

sin y , K/cVp i 

Evidently, we have ^/K/tVp j < n/lfor every location p v Moreover, cf^ft) < 3 under the condition 

q°P'{-)eC\(0,l). 

Finally, note that a wide family of classical impulsive control systems (see e.g., [12]) can be 
described by the conventional controllable Euler-Lagrange or Hamilton equations (see [5]). 
Moreover, we refer to [6] for impulsive hybrid control systems and associated OCPs. Thus the 
impulsive hybrid systems of mechanical nature can also be incorporated into the modeling 
framework presented in this section. 

4. The variational approach to hybrid OCPs of mechanical nature 

An effective numerical procedure, as a rule, uses the specific structure of the problem under 
consideration. Our aim is to study the variational structure of the main OCP (5). Let 

Ti := {?(■) G C* ([£;_i,f,]) : 7 (f,_i) = c M , 7 (t t ) = a}, 

where i = 1, ..., r.. The vectors c,, where i = 1, ..., r are defined by the corresponding switching 
mechanism of a concrete hybrid system. We refer to [3,4,26] for some possible switching rules 
determined for various classes of hybrid control systems. We now present an immediate 
consequence of the classical Hamilton Principle from analytical mechanics. 

Theorem 1. Let all Lagrangians Lp,{t,q,q,u) be a strongly convex function with respect to the 
generalized variables c\ t , i = 1, ..., n. Assume that every boundary -value problem from (3) has a unique 
solution for every u(-) e Wf|Cj,(0,l). A function q u (-), where u(-) 6 UC] 0^(0,1), is a solution 
of the sequence of boundary-value problems (3) if and only if a restriction of this function on every time 
interval [f,-_if;), i = l,...,r can be found as follows 

<?"(•) = ar g mi V)er, / ' L pi {t,q(t),q(t),u(t))dt. 
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For an admissible control function u() from U we now introduce the following two 
functionals 

ft, 
T pi (q(-),z(-)) := / {L Pi (t,q{t),q(t),u(t)) - L pi (t,z(t),z{t),u(t))}dt, 

ft, 
V Pi (q(-)):= max / [L Vi {t,q{t),q{t),u{t)) - L Vt {t,z(t),z{t),u{t)))dt, 
z(-)er, ■/*,_! 

for all indexes p, £ V '. Generally, we define q u (-) as an element of the Sobolev space 
W,,' (0,1), i.e., the space of absolutely continuous functions with essentially bounded 
derivatives. Let us give a variational interpretation of the admissible solutions (?"(•) to a 
sequence of problems (3). 

Theorem 2. Let all Lagrangians L pi {t, q, q, u) be strongly convex functions with respect to the 
variables qi, i = 1, ..., n. Assume that every boundary -value problem from (3) has a unique solution for 
every u(-) e Wf|Cj„(0,l). An absolutely continuous function q" (■), where «(■) e Wf|Cj,(0,l), 
is a solution of the sequence of problems (3) if and only if a restriction of this function on 
[h—lU)t i = 1' •••' r can ^ e found as follows 

<?"(•) = ar S mm < ? (.) 6 w,V~(t,_ 1 ,f,) y P.('?(')) (6) 

Proof. Let q ll (-) G W n ' (ti—l/h) be a unique solution of a partial problem (3) on the 
corresponding time interval, where m() G Wf|Cj„(0, 1). Using the Hamilton Principle in 
every location p, inP, we obtain the following relations 

f'< 
min ^p (<?(■)) = mm lnax / [Lp {t>q{t),q{t)> u {t))~ 

Cj(-)ewT(t t -i_,ti) ' q(-)€W 1 /°(t i - 1 ,t i ) z (-)^ T i Jt i-j 

f ' L Pi {t,z(t),z(t),u(t))}dt = min / ' L p {t,q(t),q(t),u(t))dt- 

rt, ft, 

in / L„.(t,z(t),z(t),u(t))dt= L p .{t,q u (t),q u (t),u{t))dt- 

f' L pi (t,q"(t),q"(t),u(t))dt=V pi (q"(-)) 



mm 



0. 

If the condition (6) is satisfied, then q l '(-) is a solution of the sequence of the boundary-value 
problem (3). This completes the proof. □ 

Theorem 1 and Theorem 2 make it possible to express the initial OCP (5) as a multiobjective 
optimization problem over the set of admissible controls and generalized coordinates 

minimize J(q(-),u(-)) and P{q{-)) 

subject to . . 

(<?(■),"(•))£( U r,)x(wnc»(o,i)), 

i=\,...,r 
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minimize ]{q(-),u(-)) and V(q(-)) 
subject to 

(<?(■),"(■))£( U r,)x(wnc?„(o,i)), 

i=\,...,r 



(8) 



where 



and 



J : 1 



The auxiliary minimizing problems (7) and (8) are multiobjective optimization problems (see 
e.g., [16,28]). Note that the set 

rx(wncJ,(o,i) 

is a convex set. Since fo(t, -, •), t 6 [0,1] is a convex function, J(q(-),u(-)) is also convex. If 
P(-) (or V(-)) is a convex functional, then we deal with a convex multiobjective minimization 
problem (7) (or (8)). 

The variational representation of the solution of the two-point boundary-value problem (3) 
eliminates the differential equations from the consideration. The minimization problems 
(7) and (8) provide a basis for numerical algorithms to the initial OCP (5). The auxiliary 
optimization problem (7) has two objective functionals. For (7) we now introduce the 
Lagrange function [28] 

A(t,q(-),u(-), F , F3 ) := Fl J{q(-),u(-)) + m P(q(-))+ 
p 3 |/<|dist( U(=i r r,)x(wnc, 1 „(0,i)){('?(')' »(•))}/ 

where dist/y r^xfwnc 1 (01)){'} denotes the distance function 

dist (r«)x(«nq,(o,i)){(l(i"('))} := infill (?(•)/ "(•))- 

-ellc;,(o,i)xc;„(o,i)'<? G ( U r ;) x (wncl,(o,i))}. 

i=l,...,r 
We also used the following notation 

¥■•■= iFl'f'lf G R+- 
Note that the above distance function is associated with the Cartesian product 

( |J T t )x(Uf]cU0A))- 

i=\,...,r 
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Recall that a feasible point (<?*(•), «*(•)) is called weak Pareto optimal for the multiobjective 
problem (8) if there is no feasible point (q(-), «(■)) for which 

/(*(•),«(•))</(**(•),«*(•)) and ?(<?(•))< ?(«?*(•)). 

A necessary condition for (?*(•)/ M * (')) to be a weak Pareto optimal solution to (8) in the sense 
of Karush-Kuhn-Tucker (KKT) condition is that for every ^3 £ R sufficiently large there exist 
u* eTR 2 + such that 

G 9 w . )iU( . )) A(f /(? *(-),M*(-)^*^3)- (9) 

By 9(o(.) „(.)) we denote here the generalized gradient of the Lagrange function A. We refer to 
[28] for further theoretical details. If P(-) is a convex functional, then the necessary condition 
(9) is also sufficient for (q* (■), u* (■)) to be a weak Pareto optimal solution to (8). Let N be a set 
of all weak Pareto optimal solutions (q*{-), «*(■)) for problem (7). Since (q°P t (-)u°P t (-)) e N, 
the above conditions (9) are satisfied also for this optimal pair (q ? 1 (-)w P f (-)). 

It is a challenging issue to develop necessary optimality conditions for the proper Pareto 
optimal (efficient) solutions. A number of theoretical papers concerning multiobjective 
optimization are related to this type of Pareto solutions. One can find a fairly complete 
review in [20]. Note that the formulation of the necessary optimality conditions (9) involves 
the Clarke generalized gradient of the Lagrange function. On the other hand, there are 
more effective necessary conditions for optimality based on the concept of the Mordukhovich 
limiting subdifferentials [20]. The use of the above-mentioned Clarke approach is motivated 
here by the availability of the corresponding powerful software packages. 

When solving constrained optimization based on some necessary conditions for optimality 
one is often faced with a technical difficulty, namely, with the irregularity of the Lagrange 
multiplier associated with the objective functional [15,20]. Various supplementary conditions 
(constraint qualifications) have been proposed under which it is possible to assert that the 
Lagrange multiplier rule holds in "normal" form, i.e., that the first Lagrange multiplier is 
nonequal to zero. In this case we call the corresponding minimization problem regular. 
Examples of the constraint qualifications are the well known Slater (regularity) condition 
for classic convex programming and the Mangasarian-Fromovitz regularity conditions for 
general nonlinear optimization problems. We refer to [15,20] for details. In the case of a 
conventional multiobjective optimization problem the corresponding regularity conditions 
can be given in the form of so called KKT constraint qualification (see [28] for details). In 
the following, we assume that problems (7) and (8) are regular. 

Consider now the numerical aspects of the solution procedure associated with (7) and recall 
that discrete approximation techniques have been recognized as a powerful tool for solving 
optimal control problems [3,25,29]. Our aim is to use a discrete approximation of (7) and 
to obtain a finite-dimensional auxiliary optimization problem. Let N be a sufficiently large 
positive integer number and 

Si ■= {t = t i _ 1 ,t j ,...,t i =tij 
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be a (possible nonequidistant) partition of every time interval [t,_i, (,•], where i = 1, ..., r such 
that 

max \t> +1 - t>\ < Z? . 

and limjv^coff = for every /' = l,...,r. Define A,^' +1 := t{ +l - t\, ) = 0,...,N - 1 and 
consider the corresponding finite-dimensional optimization problem 

minimize } N {q N (-),u N {■)) and P N (q N {-)), 

(<? N (-),« N (-))e( U Tf)x(U N f]Cl N (0,l)), ( lt} ) 

i=\,...,r 

where J N and P N are discrete variants of the objective functionals / and P from (7). Moreover, 
T^ is a correspondingly discrete set T; and CJ nN (0, 1) is set of suitable discrete functions 
that approximate the trajectories set C„, (0, 1). Note that the initial continuous optimization 
problem can also be presented in a similar discrete manner. For example, we can introduce 
the (Euclidean) spaces of piecewise constant trajectories q ( •) and piecewise constant control 
functions u N (•). As we can see the Banach space C„ (0,1) and the Hilbert space L^([0, 1]) will 
be replaced in that case by some appropriate finite-dimensional spaces. 

The discrete optimization problem (10) approximates the infinite-dimensional optimization 
problem (7). We assume that the set of all weak Pareto optimal solution of the discrete 
problem (10) is nonempty. Moreover, similarly to the initial optimization problem (7) we 
also assume that the discrete problem (10) is regular. If P(-) is a convex functional, then the 
discrete multiobjective optimization problem (10) is also a convex problem. Analogously to 
the continuous case (7) or (8) we also can write the corresponding KKT optimality conditions 
for a finite-dimensional optimization problem over the set of variables (q N (-),it N (■)). The 
necessary optimality conditions for a discretized problem (10) reduce the finite-dimensional 
multiobjective optimization problem to a system of nonlinear equations. This problem can be 
solved by some gradient-based or Newton-type methods (see e.g., [24]). 

Finally, note that the proposed numerical approach uses the necessary optimality conditions, 
namely the KKT conditions, for the discrete variant (10) of the initial optimization problem (7). 
It is common knowledge that some necessary conditions of optimality for discrete systems, for 
example the discrete version of the classical Pontryagin Maximum Principle, are non-correct 
in the absence of some restrictive assumptions. For a constructive numerical treatment of 
the discrete optimization problem it is necessary to apply some suitable modifications of 
the conventional optimality conditions. For instance, in the case of discrete optimal control 
problems one can use so-called Approximate Maximum Principle which is specially designed 
for discrete approximations of general OCPs [21]. 

5. Mechanical example 

This section is devoted to a short numerical illustration of the proposed hybrid approach 
to mechanical systems. We deal with a practically motivated model that has the following 
structure (see Fig. 1). 

Let us firstly describe the parameters of the mechanical model under consideration: 
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iff' B 
B ' 



After switching 



Before switching 

Fig. 1. Mechanical example 

q\ it corresponds to the position of motor. 

q 2 is the position of inertia ] 2 . 

]\, ]i are the external inertias. 

]m is an inertia of motor. 

B m it corresponds to the friction of the motor. 

B\, B 2 they correspond to the frictions of the inertias ]\, ] 2 . 

k is a constant called the rate or spring constant. 

u it corresponds to the torque of motor. 

The relations for the kinetic potential energies give a rise to the corresponding Lagrange 
dynamics: 

K(t) 

1, 



1, , 1, .2 
2W + 2^<?2 



V(t) = -k( m -q 2 y 



Finally, we have 



L(q(t),q(t)) = -] m q\ 



1,1 9 

2/242- 2 fc ('?i"'?2) 



and the Euler-Lagrange equation with respect to the generalized coordinate q\ has the 
following form 

J m q\ + B m q 1 -k(q 2 (t)-q 1 (t)) = u(t) (11) 

We now considered the Euler-Lagrange equation with respect to the second generalized 
variable, namely, with respect to q 2 



d dL(q(t),q(t)) dL(q(t),q(t)) 
dt dq 2 dq 2 



-B 2 q 2 (t) 



We get the next relation 



} 2 q 2 (t) + B 2 q 2 (t) + k(q 2 (t) - q x (t)) = 
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The redefinition of the states X\ := q\, x 2 := (\\, X3 := c\i* x i '■= <?2 with X := (x\, x 2 , X3, x^) 7 
implies the compact state-space form of the resulting equation: 
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U, Xn 



(12) 



The switching structure of the system under consideration is characterized by an additional 
inertia )\ and the associated friction B\. The modified energies are given by the expressions: 
the kinetic energy: 

111 

^]mql + 2-/i<7i + 2^4i 

the potential energy: 



K(t) 



V(t) 



-k (t/i - q 2 ) 



The function of Lagrange can be evaluated as follows 

M<M) = 2^1+ 2^i + 2^~ 2 ki * qi ~ qi)2 
The resulting Euler-Lagrange equations (with respect to (]\ and to q 2 can be rewritten as 

Urn + h)iji(t) + (B m + Bi)4i(*) - k(q 2 (t) - q t (t)) = u(t) 
] 2 q 2 (t) + B 2 q 2 (t) + k(q 2 (t) - q t (t)) = 



(13) 



(14) 



Using the notation introduced above, we obtain the final state-space representation of the 
hybrid dynamics associated with the given mechanical model: 
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k =k 
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-Xi\ 
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(15) 



The considered mechanical system has a switched nature with a state-dependent switching 
signal. We put X4 = — 10 for the switching-level related to the additional inertia in the system 
(see above). 



Our aim is to find an admissible control law that minimize the value of the quadratic costs 
functional 

!(«(.)) = - J f [x r (f)QX(f) + u T (t)Ru(t)] dt — ► min (16) 



mm 
«(•) 



The resulting Linear Quadratic Regulator that has the follow form 

u°V\ t ) = -R- 1 (f)B T (f)P(f)X opt (f) 



(17) 



On Optimization Techniques for a Class of Hybrid Mechanical Systems 



159 



where P(t) is a solution of the Riccati equation (see [7] for details) 

P(t) = -(A T (t)P(t) + P(t)A(t)) + P(t)B(t)R-i(t)B T (t)P(t) - Q(t) 
with the final condition 



(18) 



P{t f ) = (19) 

Let us now present a conceptual algorithm for a concrete computation of the optimal pair 
(w°P , X°P (•)) in this mechanical example. We refer to [7, 8] for the necessary facts and the 
general mathematical tool related to the hybrid LQ-techniques. 

Algorithm 1. The conceptual algorithm used: 

(0) Select a t sw { g 0, tr , put an index j = 

(1) Solve the Riccati euqation (18) for (15) on the time intervals [0, t sw j] U t sw \, tt 

(2) solve the initial problem (12) for (17) 

(3) calculate x±(t sw i) + 10, if \ Xi(t sw i) + 10 |= e for a prescribed accuracy e > then Stop. Else, 
increase j = j + 1, inprove t sw { = t swi + At and back to (1) 

(4) Finally, solve (15) with the obtained initial conditions( the final conditions for the vector X(t sw j)) 
computed from (12) 
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Fig. 2. Components of the optimal trajectories 
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Finally, let us present the simulation results (figure 2). As we can see, the state X\ satisfies the 
switching condition x^ + 10 = 0. The computed switching time is equal to t sw { = 0.0057s. 
The dynamic behaviour on the second time-interval [0,50] is presented on the figure (2). 
The obtained trajectories of the hybrid states converges to zero. As we can see the dynamic 
behaviour of the state vector X°P t (t) generated by the optimal hybrid control u°P t ( ■ ) guarantee 
a minimal value of the quadratic functional !(•). This minimal value characterize the specific 
control design that guarantee an optimal operation (in the sense of the selected objective) of 
the hybrid dynamic system under consideration. 

6. Concluding remarks 

In this paper we propose new theoretical and computational approaches to a specific class 
of hybrid OCPs motivated by general mechanical systems. Using a variational structure 
of the nonlinear mechanical systems described by hybrid-type Euler-lagrange or Hamilton 
equations, one can formulate an auxiliary problem of multiobjective optimization. This 
problem and the corresponding theoretical and numerical techniques from multiobjective 
optimization can be effectively applied to numerical solution of the initial hybrid OCR 

The proofs of our results and the consideration of the main numerical concepts are 
realized under some differentiability conditions and convexity assumptions. These restrictive 
smoothness assumptions are motivated by the "classical" structure of the mechanical hybrid 
systems under consideration. On the other hand, the modern variational analysis proceeds 
without the above restrictive smoothness assumptions. We refer to [20,21] for theoretical 
details. Evidently, the nonsmooth variational analysis and the corresponding optimization 
techniques can be considered as a possible mathematical tool for the analysis of discontinuous 
(for example, variable structure) and impulsive (nonsmooth) hybrid mechanical systems. 

Finally, note that the theoretical approach and the conceptual numerical aspects presented in 
this paper can be extended to some constrained OCPs with additional state and /or mixed 
constraints. In this case one needs to choose a suitable discretization procedure for the 
sophisticated initial OCP and to use the corresponding necessary optimality conditions. It 
seems also be possible to apply our theoretical and computational schemes to some practically 
motivated nonlinear hybrid and switched OCPs in mechanics, for example, to optimization 
problems in robots dynamics. 
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1. Introduction 



In this chapter, we discuss the problem of systems control. This problem represents the most 
important challenge for control engineers. It has attracted the interest of several authors and 
different approaches have been proposed and tested. These approaches can all be divided 
into two categories; Linear and Nonlinear approaches. In linear approaches, the analysis and 
the synthesis are simple however the results are limited to a specified range of operation. In 
nonlinear approaches, the results are valid in a large domain however the analysis is very 
complex. We should also note that some works on feedback control are dedicated to the 
feedback linearization in order to make the models, when it is possible, linear by using a 
preliminary feedback. 

The most important and well-known methodologies about control analysis and feedback 
control are the following: PID approach, Describing function method, adaptive control, 
robust control, Lyapunov stability, singular perturbation method, Popov criterion, center 
manifold theorem and passivity analysis. 

The first step in the controller design procedure is the construction of a truth model which 
describes the dynamics of the process to be controlled. The truth model is a simulation model 
that includes the basic characteristics of the process but it is too complicated to be used in 
the control design. Thus, we need to develop a simplified model to be used instead. Such a 
model is defined by Friedland (Friedland, 1991) as the design model. The design model 
should capture the essential features of the process. 

In order to describe the behavior of the process, a continuous dynamic system constituted 
by a finite set of ordinary differential equations of the following form is used: 

x=F[t,x(t),u(t)] X(t ) = x 
y(t) = H[t,x(t),u(t)] 

where the state x e R", the input u e R m , the output y e Rf, and F and H are vector-valued 
functions with F : RxR" xR» -> R" and H : RxR" *R m -> RP. 
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A second kind of used model is the discrete dynamic system defined by a finite set of 
difference equations: 

x{k + l) = F[k,x(k),u(k)] *(k )=x 
y{k) = H[k,x{k),u{k)] 

where x(k) = x(kh), u(k) = u(kh), h is the sampling time interval, and k > is an integer. 

The objective of this chapter is to propose a new strategy for control design using 
optimization method which is suitable for real time applications. This new methodology is 
based on neural network which is the classical approach to treat practical results using 
experimental tests. In order to illustrate this methodology and its applications, we will 
present an example of the intake air manifold control in a Diesel internal combustion engine. 

The chapter is divided as follows: In the second section a short overview of classical control 
methods is presented. In the third section a new methodology for control is proposed. In the 
fourth section, we present the application of the new control methodology to the Diesel 
engine. And finally, we end this chapter with our conclusions and remarks. 

2. Overview of classical control methods 

A main goal of the feedback control system is to guarantee the stability of the closed-loop 
behavior. For linear systems, this can be obtained by adapting the control parameters of the 
transfer function which describes the system in a way so that the real parts of its poles have 
negative values. Otherwise, Nonlinear control systems use specific theories to ensure the 
system stability and that is regardless the inner dynamic of the system. The possibility to 
realize different specifications varies according to the model considered and the control 
strategy chosen. Hereafter we present a summary of some techniques that can be used: 

2.1 Theory of Lyapunov 

Lyapunov theory is usually used to determine the stability properties at an equilibrium 
point without the need to resolve the state equations of the nonlinear system. Let us 
consider the autonomous non-linear system 

x = F(x) (3) 

where x e R" is the state variable of the system and F is a smooth vector field. Assume that 
there is a function V defined as follows: 

V : R" — > R + so that V(x) = •» x = and lim V(x) = +co If the derivative of V along the 
trajectories of (3) is so that : 

V=^VV(x),F(x)>- <0 for all x*0 (4) 

where V designates the gradient and -<.,.>- denotes the scalar product, than the system (3) is 
globally asymptotically stable. This is the Theorem of Lyapunov (Hahn, 1967). This 
approach has been improved in the principle of Krosoviskii-LaSalle (Hahn, 1967). In fact, it 
is shown that the condition given by (4) can be relaxed to 
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V =■< VV(x),F(x) >-< for all x (5) 

under the hypothesis that the more invariant set, by (3), included in 

n = {xeR7V = 0} (6) 

is reduced to the origin. 

These two theorems are the base of a large number of results on analysis of stability for 
nonlinear systems. In fact, the theory of Lyapunov- Krosoviskii-LaSalle is fundamental and 
is the base of this analysis. In the literature, this theory can have various versions according 
to the nature of the problem, for instance, for discrete models, stochastic systems or partial 
differential equations. 

In addition to the methodologies developed before, the theory is used to describe the control 
problems. The use of this theory is illustrated by the following result of feedback stabilization. 

Let us consider the following controlled system 

x = F(x) + u-G(x) (7) 

where x e R" is the state, ueR is the control variable, F and G are smooth vector fields. 
Assume there is V a Lyapunov function so that 

V=<VV(x),F(x)>-<0 (8) 

Under some hypothesis is proved (Outbib, 1999) that the closed-loop system defined from 
(7) with 

u = -<VV(x) r G{x)>- (9) 

is globally asymptotically stable at the origin. A simple example to illustrate this result is the 
scalar system 

x = u (10) 

Clearly, the system verifies the hypothesis with V(x) = 1 / 2 ■ x 2 and the stabilizing 
control u = -X can be deduced. This approach has been applied to practical process (Outbib, 
2000; Dreyfus, 1962) 

2.2 Adaptive control 

The adaptive control is mainly used in cases where the control law must be continuously 
adapted due to the varying nature of the system parameters or its initial uncertainties. 

Let us consider the following non linear system 

x = F(x,0) (11) 

Where x denotes the state variable of the system and designates a parameter. The 
adaptive control is used in the situation where the parameter 6 is not known or can change. 
For example, let us consider the scalar classical system: 
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X = ■ x 2 + u 



(12) 



If is known the system (12) can be globally asymptotically stable using a control law of the 
form u = -0 ■ x - k(x) , where k is any smooth scalar function defined as follow: k(x)x > for 
X*Q . 

The certainty-equivalent controller is defined by 



-0-x 2 -k(x) 
= w 



(13) 



where w is the update law. 

Let V be the Lyapunov function defined by: 

V(x,0) = ^x 2 + ^-(0-8) 2 (14) 

with a > . The derivative of the closed-loop system defined from (12) and (13): 

\u = -{6-6)-x 2 -k(x) 

h ' (15) 

| = w 

is given by 

V{x,6) = -x-k{x)-{ y 0-6)-x i +a-{ y 6-0)-w (16) 

Now, if we let w = (1 / a) • x , we get 

V(x,0) = -x-k(x)<O for allx fl „ 

This implies that (x,0) is bounded and x converges to zero and ensures that the system 
(12) can be stabilized at the origin. 

2.3 Sliding mode control 

The Russian school developed the methodology of sliding mode control in the 1950s. Since 
this time, the technique has been improved by several authors (Slotine, 1984; Utkin, 1992; 
Sira-Ramirez, 1987; Bartoloni, 1989; Outbib & Zasadzinski, 2009). This approach is 
applicable to various domains, including aerospace, robotics, chemical processes, etc. 

The sliding mode control is a variable structure control method. Its principle is to force the 
system to reach and to stay confined over specific surfaces where the stability can be 
ensured, and that is based on discontinuous control signal. 

In order to illustrate the approach based on variable structure control, we now present a 
simple example. Let us consider the scalar system defined by: 
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x = u (18 ) 

Our goal is to propose a control law of the form u = u(x) so that lim x(t) = lim x(t) = . 

X— >+GO X— >+0O 

Clearly, the system (18) can be globally asymptotically stable using a control law of the form 
u = f(x,x) . In fact, one can use for instance u = -x - x . 

Now a simple analysis can show that there is no linear control law, of the form u = ax , 
which makes the system globally asymptotically stable at the origin. But, if we consider a 
state feedback that commutes between two linear laws of the form: 

■ x if x ■ x > 

■ x if x ■ x < 

than the system can be globally asymptotically stable using appropriate values for a x 
and a 2 . 

2.4 Optimal control 

The objective of the optimal control method is to search for the best dynamic course which is 
capable of transporting the system from an initial state to a final desired state at minimum 
cost. An example of its various applications can be found in the satellite control. More 
precisely, the optimal control technique can be defined as follows: 



Let us consider the following system: 

x = F(t,x(t),u(t)) 



(20) 



where x s R" designates the state variable and u e R m is the control variable. 
f :RxR n x R m — > R" is a smooth vector-valued function .The optimal control is to find a 
suitable dynamic control u(t) which allows the system to follow an optimal trajectory x(t) 
that minimizes the cost function : 

] = £H(t,x(t)At)) (21) 

Several approaches have been used to resolve this problem. Among these approaches we 
can cite the variational calculus (Dreyfus, 1962), the maximum principle of Pontryagin 
(Pontryagin, 1962) or the procedure of dynamic programming method of Bellman (Bellman, 

1957). 

Let us consider a simple example such as the hanging pendulum. The equation describing 
the movement of the pendulum angular position under an applied torque a is given by: 

0(t) + A-0(t) + m 2 -0(t) = a(t) 
0(0) = 1 0(0) = 2 

where designates the angular position at time t . Let x be the system state variable 
x = [0(t),0(t)~\ , we can write : 
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■v(') = l . X \ \ = f{x,a) (23) 

x 2 -a> ■x 1 +a) 

Therefore the optimal control goal can be to minimize the time interval r , in order to reach 
the state values x(r) = . 

2.5 Robust control 

The objective of robust control is to find a control law in such a way that the response of the 
system and the error signals are maintained to desired values despite the effects of 
uncertainties on the system. The uncertainties sources can be any disturbance signals, the 
measurement noise or the modeling errors due to none considered nonlinearities and time- 
varying parameters. 

The theory of robust control began in the 1970s and 1980s (Doyle, 1979; Zames, 1981) with 
some aircraft applications. Actually, its applications concern different domains (aerospace, 
economics, ...). 

3. New algorithm for Optimized Nonlinear Control (ONC) 

The objective of this methodology is to propose a system optimized dynamic control which 
can be used in real time control applications. The proposed methodology (Omran, 2008b) 
can be divided into five steps: 1) Modeling process, 2) Model validation, 3) Dynamic 
optimization process, 4) Creation of a large database of the optimal control variables using 
the dynamic optimization process, 5) The neural network controller. 

In the next sub-sections, we present the different methodology steps and we explain its 
application using the example of the Diesel engine system. 

3.1 Modeling process 

The general equations which describe the functioning of a system can be expressed using 
the following form (Rivals, 1995): 

X = F(X,I,u,t) 

(24) 
Y = g(X,I,u,t) K ' 

Where F and g are nonlinear functions, X is the system state variables, I is the inputs 
variables, u is the control variables to be tuned and Y is the output. 

3.2 Experimental validation 

In this phase we used specified experimental data to identify the model parameters 
used in the modeling process (models of representation: transfer function or neural 
networks, models of knowledge,...), and than we used dynamic experimental data to test 
the model responses accuracy and its validation. This step is classic in any modeling 
process. 
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3.3 Offline dynamic optimization 

In this step we present the optimization technique of the control variables over dynamic 
courses and we define the objective function to be used. The question that we should ask is 
the following: In response to dynamic inputs I(t) which solicit the system over a certain 
interval of time [0,T], what is the optimal continuous values of the control parameters u(t) 
which minimize the cumulative production of the output Y(t). Therefore the objective 
function to be minimized can be written using the following form: 

T T 

MinJY(t)-dt = Minjg(X,I,u,t)-dt (25) 

" ; o "' o 

The optimization problem has the following equalities and inequalities constraints: 

Equalities constraints: —— = F(X,I,u,t) (26) 

at 

X 
Inequalities constraints: " (27) 

M min <u< "max 

Because the problem is nonlinear, there is no analytical solution; therefore we must 
reformulate it into its discretized form as following: 

N 

Objective function: Min V g { (X, ,!,-,«,• , t { ) (28) 

dx x —X 

Equality constraints: = F(X,I,U,t) => — ^ '- = F;(X;,I,,w,,f,) 

dt At 

X i+1 -X i -AtF(X i ,I i ,u i ,t i ) = (29) 



Inequality constraints: (30) 

M min < U i < M max 

The inequality constraints are the domain definition of the system's state and control 
variables; they are the lower and upper physical, mechanical or tolerance limits which 
assure a good functioning performance of the system and prevent the system damage. In 
our case, for example, the engine speed and the intake and exhaust pressure and 
temperature must vary between a lower and upper limit to prevent engine system damage 
or dysfunction. 

3.4 Creation of the optimal database 

The optimization problem explained previously necessitates a long computation time and 
therefore it cannot be directly resolved in real time applications, in addition the inputs 
evolution must be known beforehand which is not true in any real time applications. 
Consequently we propose to resolve the problem off line for different inputs profiles that 
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are very rich in information and variety and that cover a large area of possibility of the 
system's domain and then to regroup the found solutions (inputs profiles and optimal 
control variables) in a large database which will be exploited in the following step. 
Therefore in the created database, we will find for every input vector I(t) an output vector 
u(t) which is the optimal control variables that can be used to respond to the inputs 
solicitations. In the next section this database will be used to create a dynamic controller 
based on neural networks. 

3.5 Online neural network control 

Since a score of year, the neural networks are considered as a powerful mathematical tool to 
perform nonlinear regression. Many engineers used them as a black box model to estimate 
the system responses and they also used them in various fields of applications including 
pattern recognition, forms recognition, objects classification, filters and control systems 
(Rivals, 1995). We distinguish two main types of neural networks: feed-forward or multi- 
layers networks used for steady state processes and feedback or recurrent networks used for 
dynamic processes. We recognize to these networks the following fundamental 
characteristics (Ouladssine, 2004): They are black box models with great capacity for 
universal, flexible and parsimonious functions approximation. 

We are interested in establishing a control technique by training a recurrent neural network 
using the database created in the forth step of this methodology. The main advantage of this 
approach is essentially the capacity of developing a nonlinear controller with a small 
computation time which can be executed in real time applications. 

Between the various neural networks architectures found in the literature, the multi-layer 
perceptrons are the most popular; they are particularly exploited in system modeling, 
identification and control processes (Zweiri 2006). Many works show that the three layers 
perceptrons with one hidden layer are universal function approximation (Li, 2005); they are 
capable to approximate any nonlinear continuous function, defined from a finite multi- 
dimensions space into another, with an arbitrary fixed precision, while they require the 
identification of a limited number of parameters comparing to the development of series of 
fixed functions. In this way, they are parsimonious. 

4. Application: Optimal air control in diesel engine 

Many vehicles developers are especially interested in Diesel internal combustion engines 
because of their high efficiencies reflecting low fuel consumption. Therefore, electronics and 
common rail injection systems are largely developed and used in diesel engines along with 
variable geometry turbocharger and exhaust gas recirculation in order to reduce the 
pollution and protect the environment and the human health and to optimize the engine 
performance and fuel consumption. The future engines must respect the more restricted 
pollution legislations fixed in the European emissions standards (table 1). The particulate 
matter that are mostly emitted under transient conditions due to air insufficiency are 
expected to be reduced of a ratio 1:10 at 2010 (Euro 6) and the nitrogen oxides which are 
caused by a smaller rate of the exhaust gas recirculation due to the insufficiency of fresh air 
supplied to the engine by the compressor at low engine speed and fuel consumption 
reduction and engine performance at high speed are also supposed to be reduced to half. 
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Heavy duty vehicle 


Euro 1 
1993 


Euro 2 
1996 


Euro 3 
2000 


Euro 4 
2005 


Euro 5 
2008 


Euro 6 
2010 


Oxides nitrogen 


9 


7 


5 


3,5 


2 


1 


Carbon monoxide 


4,5 


4 


2,1 


1,5 


1,5 


1,5 


Hydro-carbons 


1,23 


1,1 


0,66 


0,46 


0,46 


0,46 


Particulate Matter 


0,4 


0,15 


0,1 


0,02 


0,02 


0,002 



Table 1. European standard of heavy duty vehicles in g/KW.h 

Actually, modern diesel engines are controlled by look up tables which are the results of a 
steady state optimization using experiments done on a test bench. Figure 1 shows a static 
chart of the fresh air flow rate that is used to control the air management system. This chart, 
as well as the entire look up tables used in the engine control, depends over two entries: the 
crankshaft angular speed and the fuel mass flow rate (Arnold, 2007). The schematic 
description of an open and closed loop control are shown in fig. 2 and 3. The inputs are the 
pedal's position X p and the engine angular speed w. The outputs are the actuators of the 
turbine variable geometry GV and the opening position of the exhaust gas recirculation 
EGR. The indication ref designates a reference value and the indication corr is its corrected 
value. P and m'c are respectively the predicted or measured intake pressure and the air mass 
flow rate entering the intake manifold. 



Air Flow rate 




Fuel flow rate 

id" ^s_ ^^-^ 

^i^— -**^ 1000 
500 Crankshaft Angular speed 

Fig. 1. Static chart of the fresh air flow rate used in the engine control schemes. 
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Fig. 2. Open loop control 
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Fig. 3. Closed loop control 

In the open loop control, the classic control of a diesel engine (Hafner, 2000) is done 
according to the diagram in fig. 2, the optimal values of the actuators are updated by 
memorized static maps. Then a predictive corrector (Hafner, 2001) is generally used in order 
to compensate the engine dynamic effects. 

In the closed loop control (fig. 3), the engine is controlled by error signals which are the 
difference between the predicted or measured air mass flow rate and the intake pressure, 
and their reference values. The controller uses memorized maps as reference, based on 
engine steady state optimization (Hafner, 2001; Bai, 2002). The influence of the dynamic 
behavior is integrated by several types of controller (PI, robust control with variable 
parameters, ...) (Jung, 2003). 

Our work proposes practical solutions to overcome and outperform the control insufficiency 
using static maps. The advantage of this approach is to be able to propose dynamic maps 
capable of predicting, "on line", the in-air cylinders filling. Therefore the optimal static maps 
in fig. 1 and 2 can be replaced by optimal dynamic ones. 

We suggest a mathematical optimization process based on the mean value engine model to 
minimize the total pollutants production and emissions over dynamic courses without 
deteriorating the engine performance. We used the opacity as a pollution criterion, this choice 
was strictly limited due to the available data, but the process is universal and it can be applied 
individually to each pollutant which has physical model or to the all assembled together. 

This optimization's procedure is difficult to be applied directly in "on line" engines' 
applications, due to the computation difficulties which are time consuming. Consequently, it 
will be used to build up a large database in order to train a neural model which will be used 
instead. Neural networks are very efficient in learning the nonlinear relations in complex 
multi-variables systems; they are accurate, easy to train and suitable for real time applications. 

All the simulations results and figures presented in this section were computed using 
Matlab development environment and toolboxes. The following section is divided to four 
subsections as follows: I Engine dynamic modeling, II Simulation and validation of the 
engine's model, III Optimization over dynamic trajectories, IV Creation of Neural network 
for "on line" controller. 



4.1 Engine dynamic modeling 

Diesel engines can be modeled in two different ways: The models of knowledge quasi-static 
(Winterbonne, 1984), draining-replenishment (Kao, 1995), semi mixed (Ouenou-Gamo, 2001; 
Younes, 1993), bond graph (Hassenfolder, 1993), and the models of representation by 
transfer functions (Younes, 1993), neural networks (Ouladssine, 2004). 
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Seen our optimization objective, the model of knowledge will be adopted in this work. The 
semi-mixed model is the simplest analytic approach to be used in an optimization process. 

The Diesel engine described here is equipped with a variable geometry turbocharger and 
water cooled heat exchanger to cool the hot air exiting the compressor, but it doesn't have 
an exhaust gas recirculation system that is mainly used to reduce the NOx emissions. 

Consequently the engine is divided to three main blocks: A. the intake air manifold, B. the 
engine block, C. the opacity (Omran, 2008a). 



4.1.1 Intake air manifold 

Considering air as an ideal gas, the state equation and the mass conservation principle gives [4]: 



dP 



(31) 



m c is the compressor air mass flow rate, rh aQ is the air mass flow rate entering the engine, 
P a , V a and T fl are respectively the pressure, the volume and the temperature of the air in the 
intake manifold and r is the mass constant of the air. m a0 is given by: 



\o=1v-m a0 ,th 



(32) 



l a0,th 



is the theoretical air mass flow rate capable of filling the entire cylinders' volume at 



the intake conditions of pressure and temperature: 

Vcyl.a.P, 



l a0,th 



4-n-r-T„ 



(33) 



V cy i is the displacement, a is the crankshaft angular speed, and r\ v is the in-air filling 
efficiency given by: 



= a n + a 1 a + a-,a 



(34) 



Where a, are constants identified from experimental data. The intake temperature T„ is 
expressed by: 



(/ ^ech)' '-c +r lech ' '-water 



(35) 



T c is the temperature of the air at the compressor's exit. T wa ter is the temperature of the 
cooling water supposed constant. r/ ech is the efficiency of the heat exchanger supposed 
constant. The temperature T c is expressed by: 



T=T n 
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(36) 



4.1.2 Engine block 

The principle of the conservation of energy applied to the crankshaft gives: 
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dfl 



dt\2 



;J(9> 



■P„-P.. 



(37) 



](6) is the moment of inertia of the engine, it is a periodic function of the crankshaft angle due 
to the repeated motion of its pistons and connecting rods, but for simplicity, in this paper, the 
inertia is considered constant. P e is the effective power produced by the combustion process: 



-r] e .m f . 



(38) 



7 



is the fuel flow rate, P c , is the lower calorific power of fuel and r\ e is the effective 



efficiency of the engine modeled by [5] 

r\ e =X- 



Ci +C-, -X+C, ■ X + c A -X-w + 



\c 5 -X ■w + c 6 ■ X-w +c 7 -X -w 



c: are constants, and X is the coefficient of air excess: 



P r is the resistant power: 



-C,a 



(39) 



(40) 



(41) 



Cr is the resistant torque. Fig. 4 represents a comparison between the effective efficiency 
model and the experimental data measured on a test bench. The model results are in good 
agreement with experimental data. 




Model at 800 RPM 
Model at 1200 RPM 
Model at 1600 RPM 
Model at 2000 RPM 
Exp. Data at 800 RPM 
Exp. Data at 1200 RPM 
Exp. Data at 1600 RPM 
Exp. Data at 2000 RPM 



10 20 30 40 50 60 70 80 90 100 

Air Excess 

Fig. 4. Comparison between the effective efficiency model results and the experimental data 
at different crankshaft angular speed. 
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4.1.3 Diesel emissions model 

The pollutants that characterize the Diesel engines are mainly the oxides of nitrogen and the 
particulate matters. In our work, we are especially interested in the emitted quality of 
smokes which is expressed by the measure of opacity (Fig. 5) (Ouenou-Gamo, 2001): 



Opacity = m 1 ■ vf" 1 ■ m" 



- . ( 1] 4 . |^ O O 



(42) 



m, are constants identified from the experimental data measured over a test bench. 




Air/Fuel Ratio 
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^5qq Crankshaft Angular 
Speed [rpm] 
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Fig. 5. Graphical representation of the opacity computed using (32) and a constant fuel flow 
rate equal to 6 g/s. 

4.1.4 System complete model 

Reassembling the different blocks' equations leads to a complete model describing the 
functioning and performance of a variable geometry turbocharged Diesel engine. The model 
is characterized by two state's variables (P„, w), two inputs ( thf , C,) and the following two 
differential equations representing the dynamic processes: 



a dt 

dfl, 2 
— — /o 
dt{2 



r-T a -(m c -m ao ) 
>h-m r P c ,-C r w 



(43) 



4.2 Model validation 

The test bench, conceived and used for the experimental study, involves: a 6 cylinders 
turbocharged Diesel engine and a brake controlled by the current of Foucault. Engine's 
characteristics are reported in table 2. 
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Stroke [mm] 


145 


Displacement [cm3] 


9839.5 


Volumetric ratio 


17/1 


Bore [mm] 


120 


Maximum Power [KW] 

at crankshaft angular speed [rpm] 


260 
2400 


Maximum torque [daN.m] 

at crankshaft angular speed [rpm] 


158 
1200 


Relative pressure of overfeeding [bar] 


2 



Table 2. Engine Characteristics 

Different systems are used to collect and analyze the experimental data in transient phase 
and in real time functioning: - Devices for calculating means and instantaneous measures, - 
a HC analyzer by flame ionization, - a Bosch smoke detector and - an acquisition device for 
signal sampling. The use of these devices improves significantly the quality of the static 
measures by integration over a high number of points. 

Fig. 6 and 7 show a comparison between two simulations results of the engine complete 
model and the experimental data. The inputs of the model are the fuel mass flow rate 
and the resistant torque profiles. The output variables are: the pressure of the intake 
manifold P„ the crankshaft angular speed a and the opacity characterizing the engine 
pollution. The differential equations described in section 4.1.4 are computed simultaneously 
using the Runge-Kutta method. The simulations are in good agreement with the 
experimental data. 
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Fig. 6. Simulation 1: Comparison between the complete engine model and the experimental 
data measured on the test bench. 
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Fig. 7. Simulation 2: Comparison between the complete engine model and the experimental 
data measured on the test bench. 

4.3 Optimization process 

4.3.1 Problem description 

When conceiving an engine, engines developers have always to confront and solve the 
contradictory tasks of producing maximum power (or minimum fuel consumption) while 
respecting several pollution's constraints (European emissions standard). We are only 
interesting in reducing the pollutants emissions at the engine level, by applying the optimal 
"in-air cylinders filling". Consequently, the problem can now be defined; it consists in the 
following objective multi-criteria function: 



Maximize "Power" 
Minimize "Pollutants" 



(44) 



This multi-objective optimization problem can be replaced by a single, non dimensional, 
mathematical function regrouping the two previous criteria: 



/=-f^--*+zj 



POU: 



Poll 



-it 



(45) 



P is the engine effective power, Poll is a type of pollutant, and the indication max 
characterizes the maximum value that a variable can reach. The integral represents the heap 
of the pollutants and power over a given dynamic trajectory. This trajectory can be, as an 
example, a part of the New European Driving Cycle (NEDC). 

In this chapter we will only use the opacity as an indication of pollution seen the simplicity 
of the model and the priority given to the presentation of the method, but we should note 
that the optimization process is universal and it can involve as many pollution's criteria as 
we want. The function "objective" becomes: 
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f = -\^ dt + \ 



Op 

°Pma: 



■dt 



(46) 



4.3.2 Formulation of the problem 

The problem consists therefore in minimizing the following function "objective" over a 
definite working interval [0, t] : 



■\le-m f -dt 



max 

m. 



Op i " ' 



(47) 



Under the equalities constraints representing the differential equations of the engine block 
and the intake manifold: 



dP 
V ■ — — = t • T ■ ( tH — th i 

v a ,, ' ± a V c " l aol 



— -fc> 
dt{2 



(48) 



: V e -m r P a -C r -w 



And the inequalities constraints derived from the physical and mechanical limits of the air 
excess ratio, the intake pressure and the crankshaft angular speed: 



X is given by: 



15 < X < 80 

9.5.10 4 <P a <30.10 4 [Pa] 

83<a<230 [rd / s] 



(49) 



I a + fljO + a 2 o I ■ Ncyl. Vcyl.o.P a 



4-ii-m 



(50) 



/ 



The variables of the optimization's problem are w, P a and m' c , and the inputs are C r and m'f. 

We should note that we intentionally eliminate the exhaust and turbo-compressor blocks 
from the equalities constraints because we are interesting in obtaining the optimal "in-air 
cylinders filling" m' c without being limited to any equipments such as the variable geometry 
turbo-compressor early described. In other words, we can consider that we have replaced 
the turbo-compressor by a special instrument that can deliver to the intake manifold any 
value of the air mass flow rate that we choose and at any time. Later, in the conclusion, the 
devices that can provide these optimal values are briefly discussed. Also we should note 
that we will use the complete engine model of the existing turbocharged diesel engine as a 
comparison tool, to prove the validity of our proposed optimal control and the gain in the 
opacity reduction. 
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4.3.3 Problem discretization 

There is no analytic solution to the problem previously formulated; therefore there is a 
necessity to reformulate it in its discretized form. The integrals in the function "objective" 
become a simple sum of the functions computed at different instant tf. 



f=Zf.=fl + f2+--fN 



(51) 



N is the number of the discretized points, it is the size of the unknown vectors a> , P a and m c . 
h is the step of discretization. Using the Taylor development truncated at the first 
differential order, the equalities constraints become: 



h 



1 a(i+l) 
,2 



"(•) y 



{ 



And the inequalities constraints: 






15 < /L (l) < 80 

9.5.10 4 <P a{i) <30.10 4 [Pa] 

83 < a> (l) < 230 [rd / s] 



(52) 



(53) 



4.3.4 Solution of the optimization problem 

The optimization problem under equality and inequality constraints can be described using 
the following mathematical form: 

Min{f(X)} 

X = (x 1 ,x 2 ,...x n ) 

Under Constraints (54) 

h,(X) = i = l,...,m 

gi(X)<0 i = l P 

Where f(X) is the objective function, h(X) the equality constraints and g(X) the inequality 
constraints. The easiest way to resolve this problem is to reduce it to a problem without 
constraints by creating a global objective function €>(X, h(X), g(X)) which regroups the 
original objective function and the equality and inequality constraints (Minoux, 1983). 

Therefore we will use a global objective function that regroups: The function "objective", the 
equalities constraints with Lagrange multipliers, and the inequalities constraints with a 
penalty function. The final objective function becomes (Minoux, 1983): 



L (X, X) = f(X) + £ A, * ht(X) + r.£fe,(X)J 2 



(55) 
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r = r , k is the number of iteration that must tend toward the infinity, and ro = 3. The 
problem will have m additional unknown variables (Lagrange's multipliers A;) to be 
determined along with the engine's variables. The algorithm of the minimization process 
adopted here is the Broyden-Fletcher-Goldfarb-Shanno B.F.G.S. that sums up as follows: 



1. 

2. 



To start by an initial solution X^ ' . 

To estimate the solution at the k iteration by: X k+1 = X k - a k ■ D k ■ Vf ( X k ) , with X is a 

vector regrouping the optimization variables, a k is a relaxation factor, \D k ■ V/(X i .)J 

represents the decreasing direction of the function, D k is an approximation of the 

Hessian matrix. 

To verify if the gradient's module of the objective function at the new vector X is under 

a certain desired value (» lO 2 ). If it is true then this solution is the optimal solution, end 

of search. Otherwise increment k and return to the stage 2. 



4.3.5 Results and discussion 

We applied the optimization process explained in the previous section to two different 
profiles of inputs variables (fuel mass flow rate and resistant torque). The time step of 
discretization h is equal to 0.01s and the time interval is equal to 3 sec, each problem has 900 
unknown variables {w, P„ andm c } with 598 equalities constraints and 1800 inequalities 
constraints. Fig. 8 and 9 show a comparison between the results of the optimization process 
and the simulation's results of the engine's complete system model for the same input 
values and at fixed position of the turbine variable geometry (completely open, GV = 0). The 
optimization's results show that we need significantly more air mass flow rate entering the 
intake manifold and higher intake pressure to reduce the opacity, while the real 
turbocharged diesel engine is not capable of fulfilling these tasks. 



Engine Complete Model's simulation 

Optimization process' results 
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Fig. 8. Comparison between the air mass flow rate and the intake pressure calculated using 
the optimization procedure and the ones simulated using the engine complete model for a 
variable fuel flow rate and a constant resistant torque equal to 1000 N.m. 
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Fig. 9. Comparison between the air mass flow rate and the intake pressure calculated using 
the optimization procedure and the ones simulated using the engine's complete model for a 
variable fuel flow rate and a variable resistant torque. 

Fig. 10 and 11 show a comparison between the simulated opacity derived from the 
optimization process and the one derived from the engine's complete system model for the 
same inputs used in fig. 8 and 9. The enormous gain in the opacity reduction proves the 
validity of the suggested optimization procedure. 



Engine Complete Model's simulation 
Optimization process' simulation 
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1.5 
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Fig. 10. Opacity reduction using the optimal values of the air mass flow rate and the intake 
pressure. Blue: simulation without optimization, red: Simulation with optimization. 
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Optimization process' simulation 
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Fig. 11. Opacity reduction using the optimal values of the air mass flow rate and the intake 
pressure. Blue: simulation without optimization, red: Simulation with optimization. 



4.4 Neural network controller 

Optimization previously done "off line", would be directly unexploited "on-line" by a 
controlling processor seen the enormous computation time that is necessary to resolve the 
optimization problem. In order to integrate the results of this optimization's procedure in a 
closed loop controller (ref fig. 3), and to be able to use it in real time engine applications, we 
suggest to use a black box model based on neurons. Neural network is a powerful tool 
capable of simulating the engine's optimal control variables with good precision and almost 
instantly. 

The neural network inputs are the fuel mass flow rate and the resistant torque, and its 
output variables are the optimal values of the air mass flow rate and the intake pressure. 
However in real time engine applications, the injected fuel flow rate is measurable, while the 
resistant torque is not. Consequently, we suggest substituting this variable by the crankshaft 
angular speed which can be easily measured and which is widely used in passenger cars 
controlling systems. 

Firstly, we need to create a large database which will be used to train the neural model, and 
which covers all the functioning area of the engine in order to have a good precision and a 
highly engine performance. The database is created using the optimization process as 
explained in subsection 4.3. 

Then we have to judicially choose the number of the inputs time sequence to be used, in 
order to capture the inputs dynamic effects and accurately predict the output variables. 
With intensive simulations and by trial and error, we find out that a neural network with 
inputs the fuel mass flow rate and the crankshaft angular speed at instant (i), (i-1) and (i-2) is 
capable of precisely predicting the optimal values of the air mass flow rate and the intake 
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pressure at current instant (i). Fig. 12 describes the neural network. The network is built 
using one hidden layer and one output layer, the activation functions of the hidden layer are 
sigmoid; the ones at the output layer are linear. 




Fig. 12. The structural design of the neural network adopted in this paper for predicting the 
optimal control of the in-air filling and the intake pressure in real time applications. 

The number of neurons in the hidden layer is determined by referring to the errors 
percentage of the points which are under a certain reference value wisely chosen; the 
errors percentage (table 3) are the results of the difference between the outputs of the 
network after the training process is completed, and the desired values used in the training 
database. 

Table 3 shows the results of the neural networks with different number of neurons in their 
hidden layer, these networks are trained with the same database until a mean relative error 
equal to lO 8 is reached or maximum training time is consumed. The values in the table 
represent the percentage of the neural network results respecting the specified error 
percentage computed with respect to the reference values. 



Number of neurons of the 
hidden layer 


Error percentage 


Relative 
error 


<1 % 


<5% 


<10 

% 


110 


57.71 


88.85 


96.71 


3.6 10-5 


120 


98.428 


100 


100 


io- 8 


130 


98.734 


100 


100 


10- 8 


140 


99 


100 


100 


10- 8 



Table 3. Results of four neural networks trained using different neurons number in their 
hidden layer and the same database. 

The neural network adopted in this paper includes one hidden layer with 140 neurons and 
one output layer with 2 neurons. Fig. 13 and 14 show a comparison between the air mass 
flow rate and the intake pressure calculated using the theoretical optimization procedure, 
and the ones computed using the neural network. The results are almost identical; the mean 
relative error is IO -6 . 
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Fig. 13. & 14. Comparison between the neural network outputs and the optimal values of the 
air mass flow rate and the intake pressure. 



5. General conclusions 

We successfully developed and validated a mean value physical model that describes the 
gas states evolution and the opacity of a diesel engine with a variable geometry 
turbocharger. Then we proposed a dynamic control based on the optimal "in-air cylinders 
filling" in order to minimize the pollutants emissions while enhancing the engine 
performance. The optimization process is described in detail and the simulation results (fig. 
8-11) prove to be very promising. In addition, the control principle as described here with 
the opacity criterion can be easily applied to other pollutants which have available physical 
model. This will be the object of future publications. 



Optimized Method for Real Time Nonlinear Control 



185 



Also, in order to overcome on line computation difficulties, a real time dynamic control based 
on the neural network is suggested; therefore the optimal static maps of the fig. 2 can be 
successfully replaced by dynamic maps simulated in real time engine functioning (fig. 15). 




Fig. 15. Proposed control in closed loop 

Finally, we should note that, in this chapter, while we did find, in theory, the optimal air mass 
flow rate and intake pressure necessary to minimize the opacity, but we didn't discuss the 
mechanical equipments required to provide the optimal intake pressure and intake air flow 
rate in real time engine applications. The practical implementation of the dynamic control is an 
important question to be studied thereafter. The use of a turbo-compressor with variable 
geometry and/or with Waste-Gate, and/ or electric compressor is to be considered. 
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1. Introduction 

Almost all feedback control systems are realized using discretized (discrete-time and 
discrete-value, i.e., digital) signals. However, the analysis and design method of 
discretized /quantized (nonlinear) control systems has not been established (Desoer et al., 
1975; Elia et al., 2001; Harris et al., 1983; Kalman, 1956; Katz, 1981). This article analyzes 
the nonlinear phenomena and stability of discretized control systems in a frequency domain 1 
(Okuyama, 2006; 2007; 2008). In these studies, it is assumed that the discretization is 
executed on the input and output sides of a nonlinear element at equal spaces, and the 
sampling period is chosen of such a size suitable for the discretization in the space. Based 
on the premise, the discretized (point-to-point) nonlinear characteristic is examined from two 
viewpoints, i.e., global and local. By partitioning the discretized nonlinear characteristic into 
two sections and by defining a sectorial area over a specified threshold, the concept of the 
robust stability condition for nonlinear discrete-time systems is applied to the discretized 
(hereafter, simply wrriten as discrete) nonlinear control system in question. As a result, the 
nonlinear phenomena of discrete control systems are clarified, and the stability of discrete 
nonlinear feedback systems is elucidated. 

h 



<> 



N{e r ) 



VH 



Si 
/; 



VH 



S 2 



G(s) 



-<£ 



Fig. 1. Nonlinear sampled-data control system. 



2. Discrete nonlinear control system 

The discrete nonlinear control system to be considered here is represented by a sampled-data 
control system with two samplers, Sj, S2 and the continuous nonlinear characteristic N(-) as 



1 In the time domain analysis (e.g., Lyapunov function method), it is difficult to find a Lyapunov function 
for the discretized (severe nonlinear characteristic) feedback system. The frequency domain analysis 
will be important in cases where physical systems with uncertainty in the system-order are considered. 
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Fig. 2. Discrete nonlinear control system. 
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(a) 



Fig. 3. Discretized nonlinear characteristics. 



(b) 



shown in Fig. 1. Here, VH denotes the discretization and zero-order-hold, which are usually 
performed in A/D(D/A) conversion, and G(s) is the transfer function of the linear controlled 
system. It is assumed that the two samplers with a sampling period h operate synchronously. 

The sampled-data control system can be equivalently transformed into a discrete control 
system as shown in Fig. 2. Here, G(z) is the z-transform of G(s) together with zero-order-hold, 
and T>\ and T>i are the discretizing units on the input and output sides of the nonlinear 
element, respectively. The relationship between e and v = N c j(e) in the figure becomes a 
stepwise nonlinear characteristic on integer grid coordinates as shown in Fig. 3 (a). Here, a 
round-down discretization, which is usually executed on a computer, is applied. Therefore, 
the relationship between e + and w + is indicated by small circles (i.e. a point-to-point transition) 
on the stepwise nonlinear characteristic. Even if continuous characteristic N(-) is linear, the 
discretized characteristic V becomes nonlinear on integer grid coordinates as shown in Fig. 3 
(b) (Okuyama, 2009). 

In Fig. 2, each symbol e, u, y, ■ ■ ■ indicates the sequence e(k), u(k), y(k), • ■ ■ , (k = 0, 1,2, • • • ) in 
discrete time, but for continuous value. On the other hand, each symbol e , U , • ■ ■ indicates a 
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discrete value that can be assigned to an integer number, e.g., 

e f e {■■■ ,-37,-27,-7,0, 7, 27, 37,---}, 
« + G {■■■ ,-37,-27,-7,0, 7, 27, 37,.--}, 

where 7 is the resolution of each variable. Here, it is assumed that the input and output signals 
of the nonlinear characteristic have the same resolution in the discretization. In the figure, 
e and u also represent the sequence e (k) and u (k). Without loss of generality, hereafter, 
7 = 1.0 is assumed. Thus, the input and output variables of the nonlinear element can be 
considered in the set of integer numbers, i.e., 

e + (fc), u + (fc) G Z 

Z = {■■■ ,-3,-2,-1,0,1,2,3,- ■■}. 

3. Equivalent discrete-time system 

In this study, the stepwise and point-to-point nonlinear characteristic is partitioned into the 
following two sections: 

N d (e)=K(e + v{e)), < K < 00, 

\v{e)\ < v < 00, (1) 

for \e\ < £, and 

N d (e) = K{e + n(e)), < K < 00, 

\n(e)\ < Dt\e\, < a < 1, (2) 

for \e\ > £, where v(e) and n(e) are nonlinear terms relative to nominal linear gain K. Equation 
(1) represents a bounded nonlinearity which exists in a finite region. On the other hand, (2) 
represents a sectorial nonlinearity of which the equivalent linear gain exists in a limited range. 
Therefore, when we consider the robust stability "in a global sense", it is sufficient to consider 
the nonlinear term 11(e). Here, e is a threshold of the input signal e. As a matter of course, (1) 
and (2) must be satisfied with respect to the discretized value e = e because e G e. 

Based on the above consideration, the following new sequences e*^(k) and w^(k) are defined: 

e*}(k)=el l (k)+ q - Ae ^, (3) 

w„l(k)=iol(k)-* q - Ae ^. (4) 

where q is a non-negative number, e m (k) and 10 m (k) are neutral points of sequences e (k) and 

w\k), 
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Fig. 4. Nonlinear subsystem. 
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Fig. 5. Equivalent feedback system. 

and Ac + (fc) is the backward difference of sequence e^(k), 

Ae + (fc) = e + (fc)-e + (/c-l). 



(7) 



The relationship between equations (3) and (4) in regard to the continuous values is shown by 
the block diagram in Fig. 4. In this figure, S is defined as 



8(z) 



2 1-z- 1 
h' 1 + z- 1 ' 



(8) 



Equation (8) corresponds to the bilinear transformation between z and S. Thus, the loop 
transfer function from w* to e* can be given by F(oc, q, z), as shown in Fig. 5, where 



F(oc,q,z) 



(l + qS(z))KG{z) 
1 + (1 + aqS(z))KG(z)' 



(9) 



and r' , d' are transformed exogenous inputs. Here, the variables such as w* , u' and y' written 
in Fig. 5 indicate the z-transformed ones. 

In this study, the following assumption is provided on the basis of the relatively fast sampling 
and the slow response of the controlled system. 

[Assumption] The absolute value of the backward difference of sequence e(k) is not more 
than 7, i.e., 

|Ae(fc)| = |e(*)-c(*-l)|<7. (10) 

If the condition (10) is satisfied, Ae + (fc) defined by (7) is exactly ±7 or because of the 
discretization. That is, the absolute value of the backward difference can be given as 

|Ae + (fc)| = \e f (k)-e f {k-l)\ =7 or 0. □ 

This assumption will be satisfied in the following examples. 
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4. Norm conditions 

In this section, some lemmas for norm conditions are presented. Here, in regard to (2), the 
following new nonlinear function is defined. 

/(e) := n(e) + oce. (11) 

When considering the discretized output of the nonlinear term, iv* = n(e + ), the following 
expression can be given: 

f{e f (k))=w f (k) + ote f (k). (12) 

From inequality (2), it can be seen that the function (12) belongs to the first and third 
quadrants. Considering the equivalent linear characteristic, the following inequality can be 
defined: 

< m ■■= f -^y- < 2«- (13) 

When this type of nonlinearity /3(/c) is used, inequality (2) can be expressed as 

w f {k) = n(e + (fc)) = (j8(fc) - a)e f (k). (14) 

For the neutral points of e (k) and w (k), the following expression is given from (12): 

\(f(e*(k)) +f(e f (k- 1))) = zoi(k) + «+(*)• (15) 

Moreover, equation (14) is rewritten as 

«&(*) = GS(fc) - «)£(*)• 

Since | e„ (k) \ < \e m (k)\, the following inequality is satisfied when a round-down discretization 
is executed: 

\wt(k)\ < *\et(k)\ < a\e„,(k)\. (16) 

Based on the above premise, the following norm inequalities are examined (Okuyama et al., 
1999; Okuyama, 2006). 

[Lemma-1] The following inequality holds for a positive integer p: 

Wm{ k )h,V < a \\ e m( k )hp < a|km(fc)|| 2 ,p- (17) 

Here, || • ||2,» denotes the Euclidean norm, which can be defined by 

/ v X 1/2 

II^WI|2, P := E-t 2 (fc) 
\k=l 

(Proof) The proof is clear from inequality (16). □ 

[Lemma-2] If the following inequality is satisfied in regard to the inner product of the neutral 
points of (12) and the backward difference (7): 

{w f m (k) + ae f m (k),Ae f (k)) p >0, (18) 



192 



Applications of Nonlinear Control 



the following inequality can be obtained: 

II *+ / 1 \ II ^ II *"t /T \ II 

IK, !W\\2,p < a\\ e mW\\2,p 
for any q > 0. Here, { ■ , ■ ) p denotes the inner product, which can be defined as 



(19) 



(x 1 (k),x 2 (k)) p =^x 1 (k)x 2 (k). 
k=l 

(Proof) The following equation is obtained from (3) and (4): 



"ic«iii P -iK T mii2, P = ^ 






2./' 



+ Ae + (fc) 

w m {k)-ocq — - — 



2,/' 



(20) 



Thus, (19) is satisfied by using the left inequality of (17). □ 

In regard to the input of n*(-), the following inequality can be obtained from (20) and the 
second inequality of (17) as follows: 



\Wm{k)\\ 2 , P < ct\\e* m (k)\\ 2 , v , 



when inequality (18) is satisfied. 



(21) 



5. Sum of trapezoidal areas 

The left side of inequality (18) can be expressed as a sum of trapezoidal areas. 
[Lemma-3] For any step p, the following equation is satisfied: 

a(p) := (<(fc) + <(fc),A e + (fc)) p = \ £(f(e f (k))+f(eHk-l)))Ae*(k). 



(22) 



k=\ 



(Proof) The proof is clear from (15). □ 



In general, the sum of trapezoidal areas holds the following property. 

[Lemma-4] If inequality (10) is satisfied in regard to the discretization of the control system, 

the sum of trapezoidal areas becomes non-negative for any p, that is, 

a(p) > 0. (23) 

(Proof) Since f(e*(k)) belongs to the first and third quadrants, the area of each trapezoid 



T(*):= hf(e f (k))+f(e\k-l)))Ae f (k) 



(24) 



is non-negative when e(k) increases (decreases) in the first (third) quadrant. On the other 
hand, the trapezoidal area r(fc) is non-positive when e(k) decreases (increases) in the first 
(third) quadrant. 

Strictly speaking, when (e(k) > and Ae(fc) > 0) or (e{k) < and Ae(fc) < 0), r(fc) is 
non-negative for any k. On the other hand, when (e(k) > and Ae(fc) < 0) or (e(k) < 
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and Ae(fc) > 0), r(k) is non-positive for any k. Here, Ae(k) > corresponds to Ae f (fc) = 
7 or (and Ae(fc) < corresponds to Ae (k) = — 7 or 0) for the discretized signal, when 
inequality (10) is satisfied. 

The sum of trapezoidal area is given from (22) as: 

a{p) = £ T(k). (25) 

fc=l 

Therefore, the following result is derived based on the above. The sum of trapezoidal areas 
becomes non-negative, c(p) > 0, regardless of whether e(k) (and e + (k)) increases or decreases. 
Since the discretized output traces the same points on the stepwise nonlinear characteristic, 
the sum of trapezoidal areas is canceled when e(k) (and e (k) decreases (increases) from a 
certain point (e + (fc), _f(e + (k) ) ) in the first (third) quadrant. (Here, without loss of generality, the 
response of discretized point (e + (fc),/(e + (ft:))) is assumed to commence at the origin.) Thus, 
the proof is concluded. □ 

6. Stability in a global sense 

By applying a small gain theorem to the loop transfer characteristic (9), the following robust 
stability condition of the discrete nonlinear control system can be derived. 

[Theorem] If there exists a q > in which the sector parameter a. in regard to nonlinear term 
n ( • ) satisfies the following inequality, the discrete-time control system with sector nonlinearity 
(2) is robust stable in an £ 2 sense: 

.<,K^):=^ QV+ ^ 2n2V2+ j L ;y 2){(1 + LJ)2 + V2} , (26) 

Vo; G [0, CO c ], co c '■ cutoff frequency 

when the linearized system with nominal gain K is stable. Here, fi(o>) is the distorted 
frequency of angular frequency (V and is given by 

iwhs —, ^ -2. f wh \ 



<He^)=/fiM= 7 -tan^ — j, /=</=! (27) 

In addition, U(to) and V(co) are the real and the imaginary parts of KG(& ), respectively. 

(Proof) Based on the loop characteristic in Fig. 5, the following inequality can be given in 
regard to 2 = e! wh : 

Il4(z)ll2,p< ci||r; H (2)|| 2 ,p + c 2 ||< ! (z)|| 2/P + sup|F(a,^,z)|- \\w^{z)\\2, v - 

2=1 

Here, r m (z) and d' m (z) denote the z-transformation for the neutral points of sequences r'(k) 
and d'(k), respectively. Moreover, C\ and c 2 are positive constants. 

By applying inequality (21), the following expression is obtained: 

l-a-sup|F(a,<j,z)| V* m (z)\i, v < Ci\\/ m (z)\\ 2 , p + c 2 \\d' m (z)\\ 2 , p . (28) 

2=1 / 
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Therefore, if the following inequality (i.e., the small gain theorem with respect to ^2 gains) is 
valid, 

\F(cx,q,e' a ' h )\<l/a, (29) 

the sequences e* n (k), e m (k), e(k) and y(k) in the feedback system are restricted in finite values 
when exogenous inputs r(k), d(k) are finite and p — > 00. 

By substituting (9) into inequality (29), the following is obtained: 



(l+jqd{Lo))KG(e' u ' h ) 
f (1 + jtxqCl(<v))KG{eJ ' h } 



< -. (30) 



From the square of both sides of inequality (30), 

a 2 {l + q 2 d 2 ){U 2 + V 2 ) < (1 + U - aqClV) 2 + (V + cxqClU) 2 

Then, 

a 2 {U 2 + V 2 ) + IctqCLV - {{1 + U 2 ) + V 2 } < 0. (31) 

Consequently, as a solution of inequality (31), 

-qCtV + y>q 2 n 2 V 2 + {U 2 + V 2 ){{1 + U) 2 + V 2 } 
a< U 2 + V 2 

can be given. □ 

Since inequality (26) in Theorem-1 is for all w (and n) considered and a certain q, the condition 
is rewritten as the following max-min problem. 

[Corollary] If the following inequality is satisfied, the discrete-time control system with sector 
nonlinearity (2) is robust stable: 

a < ii(qo,Wo) = max min t](q,co), (32) 

q co 

when the linearized system with nominal gain is stable. □ 

In this study, a non-conservative sufficient condition for the stability of discrete-time and 
discrete-value control systems is derived by applying the concept of robust stability in our 
previous paper(Okuyama et al., 2002a). The stability condition is, however, not satisfied for 
the entire area of the input of nonlinearity N(e) because of the stepwise and point-to-point 
characteristic. Even if the response seems to be asymptotic, there may remain a fluctuation 
(a sustained oscillation in the discrete time) or an offset. Of course, a divergent response 
that reaches the sustained oscillation may occur. These responses are typical nonlinear 
phenomena. The theorem (and corollary) derived here should be considered as the robust 
stability condition in a global sense. In addition, it is valid based on an assumption in the 
relationship between the sampling period and the system dynamics. However, this result will 
be useful in designing a discrete (digital, packet transmission) control system in practice. 

Naturally, the stability condition becomes that of continuous-time and continuous-value 
nonlinear control systems, when the sampling period h and the resolution 7 approach zero. 
Inequality (26) in Theorem-1 corresponds to Popov's criterion for discrete-time systems and 
contains the circle criterion for nonlinear time-varying (discrete-time) systems in a special 
case. The relationship between them will be described in the next section. 
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7. Relation to Popov's criterion 

Inequality (30) can be rewritten as follows: 



aH{a,q,ei wh ) 
\ + ecH(u,q,ei wh ) 



< 1, (33) 



where 

^ 1 + (1 - a)KG(eM) 

From (33), the following inequality is obtained: 

2a-K{H(a,,< ? ,e^ , ')} + l > 0. (34) 

Therefore, the following robust stability condition can be given: 

[ 1 + (1 + tt)KG(e/" ft ) + 2jaqn{w)KG(^ h ) ] 

\ l+(l-oc)KG(ei wh ) J ' ' 

which is equivalent to inequality (26). When a = 1 is chosen, (35) can be written as follows: 

l + »{(l+/.|n(tf))G(e^)} / (36) 

where _K m = 2J<C. In this case, the allowable sector of nonlinear characteristic N(-) is given as 

< N(e)e < K m e 2 , e £ 0. (37) 

When /; approaches zero (or a; is a low frequency), inequalities (36) and (37) are equivalent to 
an expression of Popov's criterion for continuous-time systems. 

In case of q = 0, the left side of (26) becomes the inverse the absolute value of complementary 
sensitivity function T(jcu). 

A ^ VU 2 + V* \T(ja>)\ 

On the other hand, from (35) 

' 1 + (I + cc)KG{ei wh ) 
1 + (1 - oc)KG((ei uh ) 

is obtained. Inequalities (38) and (39) correspond to the circle criterion for nonlinear 
time-varying systems. 



K { ; / ;V— rV } > (39) 



8. Validity of Aizerman's conjecture 

In the following case, Theorem-l becomes equal to the robust stability condition of the 
linear interval gain that corresponds to Aizerman's conjecture which was extended into 
discrete-time systems (Okuyama et al., 1998). 

[Theorem-2] If the right side of (32) is satisfied at the saddle point, 

m^\ = o, (40) 
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Fig. 6. Discretized nonlinear characteristic and stable sector for Example-1. 

inequality (26) of Theorem-1 becomes equal to the robust stability condition provided for a 

linear time-invariant discrete-time system. 

(Proof) This theorem can easily be proven by using the right side of (26). Then, 



dq 



-r]{q,ui)Cl{cjo)V{u]) 



vVnV + (u 2 + v 2 ){(i + u) 2 + v 2 } ' 

From (40), the following can be obtained: 

>/(<7o^o)^(woM^o) = 0. 
Obviously i](q,cv) > 0. Moreover, since < cog < n/h, CI(coq) > from (27). Then, 

V(cv ) = 
is obtained. Thus, 



V(lO> M o) 



|l + U(«o)| 



> a 



(41) 



(42) 



(43) 



(44) 



\U(co )\ 

Inequality (44) corresponds to the stability condition which was determined for the 
time-invariant discrete-time system with a linear gain, i.e., the Nyquist stability condition 
for a discrete-time system. 

Theorem-2 shows that the robust stability condition for a linear time-invariant system (the 
concept of interval set) can be applied to nonlinear discrete-time control systems, when (40) is 
satisfied. However, (32) is not always valid at the saddle point given in (40). In the following 
example, it can be shown that there are counter examples of Aizerman's conjecture extended 
into the nonlinear discrete-time systems. 



9. Numerical examples 

In order to verify the theoretical result, two numerical examples for discrete control systems 
with saturation type nonlinearity are presented. 
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(a) 



(b) 



Fig. 7. Time responses of e(k) and e (k), and phase traces {e(k), Ae(fc)) for Example-1 
(r = 1.0,2.0,3.0,4.0,5.0). 

[Example-1] Consider the following controlled system: 



G(s) 



K„ 



■6) 



s(s + l)(s + 2)' 



(45) 



where Kp = 1.0. It is assumed that the discretized nonlinear characteristic (discretized 
sigmoid, i.e., arc tangent function (Okuyama et al., 2002b) is as shown in Fig. 6. Here, the 
resolution value is chosen as 7 = 1.0. For C-language expression, it can be written as 

e + = y * (double) (int)(e/7), 
v = 0.4 * e + + 3.0 * atan(0.6 * e + ), 
v = 7 * (double) (int)(u/7), 

where (int) and (double) denote the conversion into an integral number (a round-down 
discretization) and the reconversion into a double-precision real number, respectively. 

When choosing the threshold £ = 2.0, the sectorial area of the stepwise (point-to-point) 
nonlinearity for £ < \e\ < 35.0 can be determined as [0.5,1.5] drawn by dotted lines in the 
figure. In this example, the sampling period is chosen as h = 0.1. From (26) and (32), the 
max-min value can be calculated as follows: 



max }](q, tv ) = i](q ,w ) 
1 



0.49, 



when the nominal gain K = 1.0. Hence, oc < 0.49 and the stable area is determined as 
[0.51,1.49]. Obviously, this sector contains the area bounded by the dotted lines. Thus, the 
discrete control system is stable in a global sense. 

The stability condition for linear gain K can be calculated as < K < 1.5 when the 
sampling period is h = 0.1. In this example, Aizerman's conjecture for discrete-time system is 
satisfied. Figures 7 (a) and (b) show time responses e(k), e + (fc) and phase traces (e(k),Ae(k)), 
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(a) 



(b) 



Fig. 8. Time responses of e(k), e (k), and phase traces (e{k), Ae(fc)) for Example-1 (r = 5.0, 
6.0, 7.0, 8.0, 9.0). 



Ae + (Jc) 



■ 



Ae + (it) 



(a) 



(b) 



Fig. 9. Backward difference Ae + (fc) vs. sampling period h for Example-1 (r = 3.0 and r = 9.0). 

(e (k),Ae (k)) of the discrete nonlinear control system when the reference inputs are r = 
1.0,2.0,3.0. Figures 8 (a) and (b) show those responses when r = 5.0,7.0,9.0. Although 
the responses contain sustained oscillations, they do not exceed the threshold e = 2.0. The 
input and the output of the nonlinearity lie in a parallelogram shown in Fig. 6. The robust 
stability in a global sense is guaranteed for all the reference inputs r. The above behavior can 
be estimated from the intersections of the highest gain of the sector and the stepwise nonlinear 
characteristic. Obviously, discrete-values (1.0,2.0) and ( — 1.0,-2.0) lie in the outside of the 
stable sector. 

Figures 9 (a) and (b) show the traces of backward difference Ae (k) when the sampling period 
/; increases. As is obvious from the figure, the assumption of (10) is satisfied for h < 0.12 when 
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Fig. 10. Discretized nonlinear characteristic and stable sector for Example-2. 

Ae(7c)| °~ 





(a) 



(b) 



Fig. 11. Time responses of e(k) and e (k) and phase traces {e(k), Ae(fc)) for Example-2 
(r = 1.0, 2.0, 3.0, 4.0, 5.0). 

r = 9.0, and for h < 2.0 when r = 3.0. In either case, the assumption is satisfied in regard to 
h = 0.1. 

[Example-2] Consider the following controlled system: 

K p (-s + 8)(s + 4) 



G(s) 



s(s + 0.2)(s + 16) 



(46) 



where Kp = 1.0. Here, the same nonlinear characteristic is chosen as shown in Example-1. 
When the threshold e = 1.0 is specified, the sectorial area of the stepwise nonlinearity for 
£ < \e\ < 10.0 can be determined as [0.78,2.0]. In this example, the sampling period is chosen 
as h = 0.04. The max-min value can be calculated as follows: 



max»(o,wo) = ri(qQ,cvo) = 0.45, 
1 
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(a) 



(b) 



Fig. 12. Time responses of e(k) and e + (fc) and phase traces (e(k), Ae(fc)) for Example-2 
(r = 5.0, 6.0, 7.0, 8.0, 9.0). 

when the nominal gain K = 1.4. Hence, a < 0.45 and the stable area is determined as 
[0.77,2.02]. This sector contains the above area. However, the stability region of control 
systems with linear gain K is given as < K < 6.3 when the sampling period is h = 0.04. 
Obviously, the discrete nonlinear control system corresponds to a counter example of the 
Aizerman conjecture. Figures 11 and 12 show time responses e(k), e + (fc) and phase traces 
(e(k),Ae(k)), (e + (fc), Ae + (fc)) of the discrete nonlinear control system, respectively. Although 
the nonlinear characteristic exists in the stable area for linear systems, a sustained oscillation 
is generated on account of a steep build-up characteristic in the lower side of the stable sector. 
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(b) 



Fig. 13. Backward difference Ae + (fc) vs. sampling period h for Example-2 (r = 3.0 and 
r = 9.0). 
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Figure 13 shows the traces of backward difference Ae (k) when the sampling period h 
increases. As is obvious from the figure, the assumption of (10) is satisfied for h < 0.11 when 
r = 9.0, and for h < 2.2 when r = 3.0. In either case, the assumption is satisfied in regard to 
h = 0.04. 

10. Conclusions 

This article analyzed the nonlinear phenomena and stability of discrete-time and 
discrete-value (discretized/quantized) control systems in a frequency domain. By 
partitioning the discretized nonlinear characteristic into two nonlinear sections and by 
defining a sectorial area over a specified threshold, the concept of the robust stability condition 
for nonlinear discrete-time systems was applied to the discrete nonlinear control systems. In 
consequence, the nonlinear phenomena of discrete control systems were clarified, and the 
robust stability of discrete nonlinear feedback systems was elucidated. The result described in 
this chapter will be useful in designing discrete (digital, event-driven, or packet transmission) 
control systems. 
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